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1  Overview 

The  interaction  of  laser  energy  with  organic  matrix  composites  (OMCs)  is 
important  in  a  wide  variety  of  processes  relevant  to  both  commercial  and 
defense  interests.  For  example,  lasers  are  used  in  the  machining  of  composites 
during  laser  cutting  of  bulk  panels  as  well  as  for  micromachining  detailed 
features.  Lasers  are  used  in  additive  manufacturing  processes,  both  in  fused 
deposition  modeling  and  laser  sintering,  and  similar  processes  are  desired  for 
additive  manufacturing  of  OMCs.  Finally,  high  energy  lasers  make  up  one 
component  of  emerging  threats  to  modern  weapons  systems.  Understanding 
the  fundamental  interactions  between  lasers  and  OMCs  will  lead  to  increased 
option  space  for  engineering  desired  responses. 

This  project  proposed  to  address  gaps  in  the  community’s  understand¬ 
ing  of  laser  interactions  in  OMCs  through  a  combination  of  modeling  and 
experimental  efforts.  Specifically,  the  project’s  tasks  include  the  following. 

1.  Develop  an  integrated  optical  and  continuum  model  for  thermal  and 
mechanical  properties  of  composites  under  illumination  from  high  en¬ 
ergy  laser. 

2.  Compare  modeled  and  experimental  results  as  a  function  of  interface 
nanotailoring  (modulation  of  absorption,  thermal  transport,  and  degra¬ 
dation). 
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3.  Explore  autonomous  sensing  and  responsive  concepts  based  on  nanopar¬ 
ticle  thermo-chromophores  and  smart  materials. 

A  key  challenge  in  each  of  these  tasks  is  the  multiscale  nature  of  the 
relevant  physics  as  illustrated  in  Figure  1.  Optical  energy  can  be  transmit¬ 
ted,  reflected,  or  absorbed  by  each  individual  component  of  the  composite. 
Absorbed  energy  is  efficiently  transduced  to  heat,  which  then  flows  through 
individual  components  as  well  as  across  interface  boundaries. 


Optical 
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Diffraction 


Slit  Diffraction 
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Depth 
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Figure  1:  Key  optical  and  thermal  phenomena  in  laser  interactions  with 
composites. 

As  illustrated  in  Figure  2,  bulk  average  optical  properties  are  deter¬ 
mined  by  nanoscale  features  such  as  constituent  dielectrics  and  the  nature 
of  material  interfaces.  Thus,  even  thin  coatings  on  fibers  can  have  a  sub¬ 
stantial  impact  on  the  overall  optical  properties  of  a  composite.  Similarly, 
thermal  interactions  are  important  at  the  nano-  (interfaces),  micro-  (ma¬ 
terial  properties  and  structure),  and  macro-  (bulk  transport)  scales.  The 
approach  taken  in  this  project  incorporates  both  modeling  and  experimen¬ 
tal  methods  to  understand,  quantify,  and  tailor  the  key  phenomena  which 
regulate  energy  transduction  multiscale  transport.  Importantly,  this  study 
does  not  address  nonlinear  interactions  such  as  ablation,  plasma  generation, 
temperature-dependent  optical  responses,  or  thermochemical  effects. 
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ifi 

Figure  2:  Understanding  laser  interactions  in  composite  materials  is  funda¬ 
mentally  multiscale  in  nature,  (left)  Bulk  heterogeneity  determines  overall 
optical  and  thermal  transport  in  the  material,  (center)  The  interactions  be¬ 
tween  material  properties  of  the  fiber,  matrix,  and  any  coatings  determine 
the  amount  of  optical  energy  that  is  absorbed  and  the  rate  of  thermal  trans¬ 
port,  which  is  highly  anisotropic,  (right)  The  smoothness,  relative  index  of 
refraction,  and  thermal  resistance  at  the  fiber-matrix  interface  are  nanoscale 
phenomena  that  govern  the  properties  at  larger  scales. 

2  Approach 

At  both  the  nanoscale  and  the  microscale,  computational  models  are  used 
to  help  develop  a  more  detailed  fundamental  understanding  of  the  energy 
transport  within  composites  subject  to  high  energy  laser  irradiation.  At  the 
micro-  and  macroscales,  COMSOL®  Multiphysics  ®  was  used  to  generate 
finite  element  models  (FEM)  of  simple  and  complex  model  systems  of  com¬ 
posites  under  laser  irradiation.  At  the  nanoscale  the  effects  of  both  chemistry 
and  structure  on  thermal  transport  have  been  explored  by  detailed  molecular 
dynamics  (MD)  models.  For  each  of  these  models,  experimental  techniques 
are  developed  to  provide  physical  parameters  and  validate  model  predictions. 
MD  models  track  thermal  energy  at  the  atomic  level,  and  they  are  used  to 
calculate  the  interface  thermal  conductance  across  the  fiber-matrix  interface. 
These  parameters  are  notoriously  difficult  to  measure  experimentally,  but  a 
novel  application  of  scanning  thermal  microscopy  (SThM)  was  successfully 
used  to  measure  the  relative  thermal  conductance  among  fibers  with  various 
coatings.  Similarly,  recording  the  maximum  temperature  experienced  within 
a  composite  is  important  for  validating  multiscale  models,  but  such  mea¬ 
surements  are  impossible  in  bulk  samples  using  state-of-the-art  thermometry 
techniques.  As  part  of  this  project,  a  temperature  sensor  based  on  the  ther¬ 
mal  reshaping  of  gold  nanorods  was  been  developed. 
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This  report  summarizes  the  key  hirelings  of  each  these  efforts,  identifies 
pitfalls  and  challenges,  and  outlines  potential  routes  forward.  It  is  organized 
according  to  the  key  tasks  listed  in  the  Overview.  The  appendices  comprise 
articles  published  in  peer-reviewed  journals  which  provide  detailed  results 
and  analyses  not  found  in  the  body  of  the  report. 


3  Develop  an  integrated  optical  and  contin¬ 
uum  model 

3.1  COMSOL®  Multiphysics 

To  develop  an  understanding  of  the  fundamental  constraints  on  energy  how 
within  OMCs  upon  laser  irradiation,  finite  element  analysis  (FEA)  was  em¬ 
ployed  to  analyze  optical  interactions,  opto-thermal  transduction,  and  ther¬ 
mal  transport  in  a  variety  of  model  systems.  COMSOL®  Multiphysics  is  a 
commercial  FEA  package  that  allows  for  multiple  physics  modules,  solvers, 
and  parameterizations  in  one,  two,  and  three  dimensions.  Models  were  con¬ 
structed  in  COMSOL®  in  variations  of  two  basic  geometries:  an  idealized 
2-D  system  and  several  full  3-D  systems  representing  experimentally  realized 
configurations. 

3.1.1  2-D  Idealized  Geometry 

An  idealized  system  consisting  of  a  single  cylindrical  hber  embedded  in  a 
rectangular  matrix  was  constructed  to  isolate  the  role  of  material  properties 
in  the  response  of  composites  to  laser  irradiation.  This  allows  the  study  of  the 
potential  to  tune  laser  response  by  varying  constituent  properties  in  a  com¬ 
posite.  First,  analytical  and  computational  models  for  the  optical  absorption 
by  carbon  hbers  were  used  to  determine  whether  or  not  the  directionality  of 
the  heat  input  is  significant.  These  computations  indicate  that  any  optical 
energy  which  penetrates  the  hber  surface  is  absorbed  within  tens  of  nanome¬ 
ters,  and  that  optical  energy  is  quickly  and  efficiently  thermalized  uniformly 
throughout  the  hber.  Thus,  the  2-D  multiphysics  models  are  constructed  us¬ 
ing  a  direct  heat  source  at  the  hber,  and  one  such  model  is  shown  in  Figure  3. 
The  dependence  of  resulting  temperature  profiles  on  bulk  thermal  conduc¬ 
tivities,  interface  thermal  conductance  (modeled  as  a  finite  coating),  surface 
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Figure  3:  Example  of  one  of  the  2-D  FEA  models  used  to  explore  the  depen¬ 
dence  of  energy  flow  on  constituent  properties  in  a  simple  composite  system. 


heat  transfer  coefficients  (including  convection  and  radiation),  fiber  optical 
properties,  and  input  power  were  investigated  using  this  model  geometry. 

Figure  4  shows  the  most  important  results  from  these  models,  and  Ap¬ 
pendix  C  comprises  a  conference  report  summarizing  the  most  significant 
findings,  ft  was  found  that,  for  uniform  heat  sources,  the  steady-state  tem¬ 
perature  gradient  within  the  homogeneous  matrix  does  not  depend  on  inter¬ 
face  thermal  resistance.  This  is  because  at  long  times  the  energy  balance  is 
uniquely  determined  by  the  temperature  at  the  outside  boundary;  the  rate  of 
energy  flow  out  of  the  system  depends  only  on  temperature  and  the  surface 
heat  transfer  coefficient. 

However,  two  features  of  the  temperature  profile  do  depend  on  the  inter¬ 
face  resistance.  First,  the  rate  at  which  the  resin  heats  up  for  a  uniform  heat 
source  at  the  fiber  is  slower  for  larger  values  of  the  interface  resistance.  Sec¬ 
ond,  the  maximum  stead-state  temperature  of  the  fibers  is  higher  for  larger 
values  of  the  interface  resistance.  Both  of  these  phenomena  are  a  result  of 
increased  thermal  energy  storage  in  the  fiber.  Therefore,  engineering  the  in¬ 
terface  thermal  resistance  between  the  fiber  and  matrix  allows  for  the  control 
of  heat  rates  and  temperature  distribution  between  the  fibers  and  the  matrix. 
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Figure  4:  Results  from  2-D  idealized  model  with  local  heat  source,  (top) 
Fiber  temperature  as  a  function  of  coating  thermal  conductivity  (ie,  effective 
interface  conductance)  for  several  thicknesses  of  coatings,  (middle)  Temper¬ 
ature  difference  between  the  fiber  surface  and  sample  surface  as  a  function  of 
coating  thickness  for  different  values  of  the  coating  conductance,  (bottom) 
Temperature  profile  in  the  resin  for  several  values  of  input  power. 
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3.1.2  3-D  Dogbone  Model 


While  the  idealized  model  results  give  insights  into  the  significance  of  various 
constituent  material  properties,  they  do  not  allow  for  quantitative  predic¬ 
tions  of  energy  flow  because  they  rely  on  independent  measurements  of  those 
properties.  To  connect  model  results  to  experimental  measurements,  a  full 
3-D  “dogbone”  geometry  was  created  to  analyze  an  experimentally  realized 
model  system  (see  Nanoparticle  Temperature  Sensing  below).  This  system 
consists  of  a  single  carbon  fiber  spanning  the  length  of  a  long,  thin  rod  with 
wide,  trapezoidal  tabs  at  each  end  as  shown  in  Figure  5  and  Appendix  B. 
The  temperature  profile  was  calculated  under  two  scenarios:  resistive  heating 
of  the  fiber  within  the  dogbone,  and  direct  laser  irradiation  in  cross-section. 
Resistive  heating  will  be  discussed  here,  and  laser  irradiation  will  be  left  for 
the  experiment  description  to  follow. 


Figure  5:  FEA  model  of  dogbone  sample  with  heat  source  from  a  single 
embedded  carbon  fiber. 

In  order  to  measure  specific  material  properties  for  accurate  FEA  mod¬ 
els,  the  COMSOL@model  of  the  dogbone  geometry  was  calculated  for  a 
wide  variety  of  material  properties.  The  bulk  thermal  conductivity,  effective 
surface  heat  transfer  coefficient,  input  power,  and  interface  resistance  were 
each  varied,  and  the  resulting  steady-state  temperature  profile  was  calcu¬ 
lated  in  each  case.  These  temperature  profiles  were  then  compared,  using 
an  automated  computer  algorithm  in  Mathematica,  to  experimentally  mea¬ 
sured  temperature  profiles,  and  the  parameter  set  corresponding  to  the  lowest 
square  deviation  was  chosen  as  representative  of  the  physical  system  (see  Ap- 
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pendix  B  and  the  Nanoparticlc  Temperature  Sensing  section  below).  Thus, 
the  COMSOL®  model  was  fully  parameterized  with  experimental  data  from 
the  corresponding  physical  system  without  needing  to  make  several  indepen¬ 
dent  measurements  of  each  property.  This  allows  for  the  model  to  be  used 
in  a  predictive  way  for  more  complex  geometries  (in  which  the  correspond¬ 
ing  physical  parameters  are  even  more  difficult  to  directly  determine).  This 
success  highlights  one  of  the  advantages  of  a  truly  integrated  computational 
materials  science  and  engineering  (ICMSE)  approach. 

3.1.3  3-D  Cross-sections 

Models  were  also  constructed  for  the  analysis  of  direct  heating  of  discrete 
slices  of  the  dogbone.  As  discussed  in  the  experimental  section  below,  direct 
laser  heating  at  the  end  of  an  exposed  fiber  was  used  to  circumvent  some  of 
the  challenges  involved  in  resistive  heating.  The  COMSOL@  models  helped 
give  insight  into  the  role  of  the  fiber-matrix  interface  in  the  axial  energy 
flow  under  these  conditions.  Figure  6  shows  typical  COMSOL®  results 
for  this  geometry.  As  with  the  idealized  2-D  geometry,  when  the  interface 
resistance  is  large  the  fiber  heats  up  more.  In  the  case  of  end-on  heating 
this  results  in  more  thermal  transport  along  the  axis  of  the  fiber  and  a  more 
even  distribution  of  temperature  along  the  axial  direction  in  the  matrix.  For 
low  values  of  the  interface  thermal  resistance  the  heat  is  readily  wicked  out 
of  the  fiber  and  the  temperature  profile  is  distributed  in  the  resin  near  the 
irradiated  surface. 

3.1.4  Complex  Microstructure 

Once  the  COMSOL®  Multiphysics  model  was  parameterized  and  validated 
by  detailed  temperature  measurements  in  the  dogbone  geometry,  it  was  then 
used  to  model  a  geometry  more  representative  of  structural  composites.  By 
incorporating  known  material  properties  into  the  fully  vetted  models,  predic¬ 
tive  computations  are  enabled  that  take  advantage  of  the  full  power  of  the 
multiphysics  solvers  to  obtain  a  complete  picture  of  the  energy  transport  in 
OMCs.  As  of  the  time  of  the  writing  of  this  report,  these  models  continue 
to  be  developed  for  specific  cases  of  interest.  A  general  picture  of  the  power 
of  the  fully  parameterized  model  is  illustrated  in  Figure  7.  Recent  efforts  to 
fully  quantify  the  thermal  response  to  laser  irradiation  in  OMCs  using  this 
approach  are  summarized  in  Appendix  C. 
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Figure  6:  COMSOL®  model  of  a  cross-sectional  slice  of  the  dogbone  sample 
upon  laser  radiation  of  the  fiber  at  one  end  or  on  the  side. 

3.2  Molecular  Dynamics  Modeling 

In  highly  filled,  heterogeneous  composite  materials,  heat  phenomena  are 
strongly  affected  by  barriers  to  thermal  transport  across  the  interfaces  be¬ 
tween  constituents.  A  molecular  dynamics  model  was  developed  to  investi¬ 
gate  the  thermal  transport  across  the  fiber-matrix  interface  for  a  variety  of 
carbon  fibers  in  a  cross-linked  bismaleimide  (BMI)  polymer  matrix  as  sum¬ 
marized  in  Figure  8.  A  detailed  description  of  the  methods  and  results  was 
published  in  ACS  Applied  Materials  and  Interfaces  (see  Appendix  A).  An 
ACS  Conference  Proceedings  report  is  also  included  in  Appendix  A  which 
describes  some  of  the  nuances  related  to  varying  time  scales  in  the  types  of 
laser  events  of  interest.  An  outline  of  the  MD  approach  and  summary  of  key 
conclusions  is  given  here. 

The  first  step  in  the  MD  approach  is  to  generate  the  atomic  coordinates 
for  the  fiber  and  the  resin.  For  the  BMI  resin,  monomeric  molecules  of  4,  4’- 
bismaleimido diphenyl  methane  (BMPM)  and  O,  O’-diallyl  bisphenol  A  were 
constructed  in  Materials  Studio  and  a  collection  of  these  monomers  were 
placed  at  random  to  fill  up  the  model  region.  The  monomers  were  allowed  to 
equilibrate  through  several  steps  as  outlined  in  Appendix  A,  but  the  complex 
chemical  cross-linking  was  not  modeled  via  MD. 

The  carbon  fiber  surface  was  treated  as  crystalline  graphite,  but  with 
several  modifications.  The  graphitic  orientation,  presence  of  sp-3  hybridized 
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Figure  7:  Multiphysics  modeling  of  full  microstructures,  (top)  Illustration 
of  different  physics  captured  by  the  model,  (bottom)  Model  results  for  full, 
randomized  microstructure  including  all  relevant  physics  for  interfaces  with 
high  (left)  and  low  (right)  interface  thermal  conductance. 
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Figure  8:  Interface  features  studied  by  MD. 


“defects” ,  thickness  of  the  modeled  graphite  region,  surface  functionalization, 
surface  roughness,  presence  of  carbon  nanotubes,  and  amount  of  hydrogen 
saturation  were  all  varied  in  the  models.  Each  of  these  modifications  corre¬ 
sponds  to  physical  parameters  that  may  be  tuneable  or  otherwise  amenable 
to  designing  an  autonomous  and/or  adaptive  response  in  the  presence  of 
optical  or  thermal  energy. 

In  each  case,  the  effective  thermal  conductance  across  the  interface  is  cal¬ 
culated  by  enforcing  a  fixed  heat  flux  across  the  interface  and  measuring  the 
average  temperature  as  a  function  of  position  along  the  flux.  The  interface 
thermal  conductance  is  calculated  according  to 


AAT 

where  Q  is  the  heat  flux,  A  is  the  cross-sectional  area,  and  T  is  the 
temperature. 

3.2.1  Graphitic  Orientation,  Defect  Density,  and  Hydrogen  Satu¬ 
ration 

For  unsaturated  carbon  fiber  surfaces,  the  MD  models  predict  a  thermal 
conductance  that  decreases  when  the  angle  between  the  graphite  planes  and 
the  interface  increases  (see  Appendix  A).  This  is  attributed  to  the  decrease  in 
low  frequency  modes  in  the  graphite  planes  compared  to  those  present  along 
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the  (001)  direction,  which  couple  readily  into  the  matrix.  For  saturated 
surfaces,  the  larger  atomic  number  density  provided  by  the  additional  small 
hydrogen  atoms  increases  this  conductance  for  all  orientations. 

Conversely,  the  presence  of  sp-3  defects  in  the  otherwise  sp-2  graphite 
planes  primarily  serves  to  scatter  these  low  frequency,  out-of-plane  modes 
and  has  minimal  impact  on  the  in-plane  modes.  Thus,  at  moderate  defect 
densities  (~  6  nm“3)  the  interface  conductance  of  every  orientation  becomes 
comparable  to  each  of  the  others. 

These  results  help  extend  current  knowledge  about  thermal  transport 
within  carbon  fibers  having  different  atomic  structure.  The  models  provide 
valuable  insights  into  the  influence  of  fiber  structure  on  bulk  composite  con¬ 
ductivity  through  variations  in  interface  thermal  conductance. 

3.2.2  Surface  Functionalization 

A  common  strategy  for  increasing  interactions  between  dissimilar  constituents 
in  a  composite  material  is  via  chemical  functionalization.  Within  the  MD 
system,  functionalization  is  modeled  by  enforcing  additional  constraints  dur¬ 
ing  the  model  construction  and  by  the  addition  of  energy  terms  that  reflect 
chemical  bonds  between  certain  monomers  and  the  carbon  fiber  surface.  As 
expected,  the  MD  results  show  that  the  interface  thermal  conductance  in¬ 
creases  with  increasing  surface  functionalization  density  (see  Appendix  A). 

3.2.3  Interface  Roughness  and  Carbon  Nanotubes 

As  the  roughness  of  the  fiber  surface  increases,  so  does  the  effective  area 
over  which  heat  transport  can  occur.  Thus,  in  general,  surface  roughness 
increases  the  interface  thermal  conductance.  Similarly,  the  presence  of  carbon 
nanotubes  grown  on  the  fiber  surface  extends  the  effective  region  of  carbon 
into  the  matrix  and  thereby  increases  the  thermal  conductance.  Quantifying 
this  increase  is  difficult,  however,  because  the  interface  itself  becomes  more 
ambiguous  at  the  atomic  level  (see  Appendix  A). 

3.2.4  Considerations 

There  are  two  main  gaps  or  limitations  in  the  insights  provided  by  the  MD 
simulations.  First,  the  MD  models  assume  perfect  interfaces,  even  with  vary¬ 
ing  surface  roughness  and  degree  of  functionalization.  Any  voids  at  the  fiber 
surface,  such  as  inevitably  arise  through  any  real  manufacturing  processes, 
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would  tend  to  decrease  the  average  interface  thermal  conductance.  Also, 
the  cross-linked  BMI  matrix  is  expected  to  have  higher  thermal  conductivity 
than  the  un-cross-linked  model  system.  This  should  result  in  a  higher  overall 
thermal  conductance  into  the  matrix,  and  this  would  be  especially  true  for 
low  frequency  vibrational  modes,  potentially  skewing  the  differences  arising 
from  graphitic  orientation,  functionalization,  and  bond  saturation. 

3.2.5  MD  Summary 

These  MD  results  offer  new  insights  into  the  importance  of  several  chemical 
and  structural  factors  in  determining  the  interface  thermal  conductance.  This 
conductance,  as  shown  later,  is  an  important  input  in  models  at  more  coarse 
length  scales,  and  it  can  make  a  large  difference  in  the  response  of  composite 
materials  to  laser  irradiation. 


4  Compare  modeled  and  experimental  results 

4.1  Nanoparticle  Temperature  Sensing 

A  key  requirement  for  the  successful  parameterization  and  validation  of  mul¬ 
tiphysics  models  for  composite  material  response  to  laser  irradiation  is  the 
precise  measurement  of  internal  temperatures  with  high  spatial  resolution. 
While  several  high-resolution  temperature  sensing  mechanisms  exist,  none 
of  them  allow  for  both  three-dimensional  mapping  of  temperature  and  the 
recording  of  temperature  for  ex-situ  measurement.  A  major  contribution 
of  this  project  was  to  develop  a  solution  with  those  capabilities.  To  that 
end,  the  optically  detected,  thermally  induced  shape  transformation  of  gold 
nanorods  in  response  to  heating  was  used  to  invent  a  nanoparticle-based, 
microscale  temperature  recording  technology. 

The  details  of  the  experimental  approach,  innovative  analytical  tech¬ 
niques,  and  resulting  technological  application  have  been  published  in  ACS 
Applied  Materials  and  Interfaces,  and  the  article  is  reproduced  in  Appendix  B 
for  reference,  ffowever,  several  experimental  details  have  been  omitted  from 
that  article  and  are  included  here  for  completeness. 


13 

Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 


4.1.1  Nanorod  Preprocessing 

In  calibration  samples,  test  samples,  and  dogbone  specimens  alike,  careful 
processing  of  the  nanocomposite  is  required  to  avoid  unwanted  shape  change 
in  the  nanorods.  Uncoated  gold  nanorods  begin  to  change  shape  as  soon 
as  their  temperature  exceeds  ~  10  —  15°  C.  Because  this  change  is  irre¬ 
versible,  any  exposure  to  elevated  temperatures  before  an  experiment  begins 
will  reduce  the  dynamic  range  of  potential  temperature  measurement.  While 
investigations  are  currently  under  way  to  engineer  nanomaterials  with  in¬ 
creased  temperature  capabilities,  a  preliminary  approach  to  this  problem 
was  to  identify  a  material  system  that  did  not  require  excessive  heating  in 
the  fabrication  stages. 

Several  materials  were  investigated  for  their  usefulness  as  room-temperature 
cure  systems.  Various  addition-cure  epoxies  and  silicones  were  tested  along 
with  solution-cast  thermoplastics  such  as  polystyrene  and  polyethylene.  Only 
the  epoxy  systems  enabled  facile  dispersion  of  the  CTAB-encapsulated  nanorods 
without  further  chemical  modification  beyond  the  scope  of  the  project.  Room- 
temperature  cures  of  these  epoxies  generally  result  in  a  very  brittle  material 
that  is  difficult  to  work  with.  Thus,  a  compromise  between  material  proper¬ 
ties  and  nanorod  shape  change  was  found  for  a  moderate  temperature  cure, 
and  the  resulting  cure  cycle  is  detailed  in  Appendix  B. 

4.1.2  Dogbone  Preparation 

The  preparation  of  the  dogbone  samples  is  delicate.  First,  a  tool  is  machined 
from  aluminum,  and  silicone  molds  are  then  cast  from  the  tool.  The  overall 
geometry  is  modeled  after  common  tensile  testing  specimens,  but  it  includes 
an  extra  channel  at  each  end  of  the  specimen  which  is  used  to  hold  a  single 
carbon  fiber  in  place.  Tape  tabs  are  placed  at  each  end  of  a  single  carbon 
fiber,  and  the  fiber  is  draped  across  the  mold  with  the  tabs  hanging  off  the 
end.  The  weight  of  the  tabs  provides  tension  during  the  cure.  Resin  is  mixed 
and  degassed  in  the  usual  way,  and  then  it  is  slowly  dispensed  at  each  end 
of  the  mold  using  a  wide  pipette  (ie,  “eye-dropper”). 

4.1.3  Resistive  Heating 

As  discussed  above,  it  was  demonstrated  that  direct  heating  of  the  carbon 
fiber,  such  as  through  Joule  heating,  is  a  suitable  surrogate  for  laser  heating 
when  the  matrix  is  transparent.  Thus,  many  of  the  experiments  used  for 
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comparison  with  the  FEA  models  were  heated  by  passing  a  current  through 
the  fiber  embedded  in  the  dogbone.  This  is  accomplished  by  placing  a  small 
amount  of  silver  paint  directly  into  the  silicone  mold.  The  paint  is  dabbed 
onto  the  ends  of  the  fiber  before  the  resin  is  added.  After  curing,  the  spot  of 
silver  paint  forms  a  contact  pad  at  either  end  of  the  specimen.  More  silver 
paint  is  used  to  connect  fine  wires  to  the  ends  of  the  sample,  and  these  wires 
are  then  used  to  inject  current  directly  into  the  fiber. 

Often  the  electrical  resistance  of  these  contacts  is  similar  to,  or  even 
greater  than,  the  resistance  through  the  fiber.  For  the  2  inch  long  (PAN) 
carbon  fibers  used  in  these  experiments  the  total  resistance  of  the  circuit 
was  typically  10  —  30  kO.  One  consistent  difficulty  in  resistively  heating  the 
dogbone  samples  arises  from  the  GTE  mismatch  between  the  carbon  fiber 
and  the  host.  In  many  samples  (>  50%)  the  thermal  strain  in  the  host  caused 
the  fiber  to  break,  which  breaks  the  electrical  circuit  and  prematurely  arrests 
the  heating.  The  results  published  in  Appendix  B  were  those  lucky  few  that 
survived  the  resistive  heating,  and  so  may  be  biased  toward  those  with  poor 
mechanical  interfaces  between  the  fiber  and  resin. 

4.1.4  Direct  Laser  Heating 

To  avoid  the  complications  arising  from  resistive  heating  of  the  dogbones,  an 
experiment  was  designed  to  realize  direct  laser  irradiation  of  the  fiber  within 
the  matrix.  A  diode  laser  with  A  =  1.1/mi  was  focused  to  a  nearly  diffraction 
limited  spot  size  (~  2/irn  diameter  gaussian  beam).  The  beam  was  aligned 
to  the  fiber  using  an  IR  imaging  camera  in  two  ways  as  shown  in  Figure  9. 

First,  to  create  a  heat  source  directly  within  the  fiber,  the  laser  was 
focused  onto  the  exposed  end  of  the  fiber  and  the  temperature  profile  was 
mapped  using  the  gold  nanorod  sensors.  As  predicted  by  the  COMSOL® 
models,  different  interface  coatings  gave  rise  to  different  temperature  profiles 
in  this  configuration  due  to  the  differences  in  axial  heat  transport.  Figure  10 
shows  two  nanorod  temperature  maps:  one  with  no  coating,  and  one  with 
an  insulating  (AI2O3)  coating. 

Second,  the  samples  were  polished  on  one  side  until  the  fiber  surface 
was  nearly  even  with  the  matrix  surface.  A  very  thin  layer  (<  0.5  mm)  of 
matrix  material  was  left  between  the  fiber  and  sample  surface.  This  was 
done  to  avoid  any  direct  heating  of  the  matrix  by  the  laser  while  maintain¬ 
ing  the  optical  properties  of  the  fiber-matrix  interface.  Unfortunately,  these 
measurements  had  to  be  suspended  due  to  lack  of  available  equipment  and 
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Figure  9:  Laser  heating  geometries  for  decoupling  interface  effects,  (left) 
End-on  irradiation,  (right)  Edge  irradiation. 
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Figure  10:  Spectroscopically  measured  temperature  maps  for  end-on  laser 
heating  of  a  carbon  fiber  with  (left)  and  without  (right)  an  insulating  coating 
of  AI2O3.  All  spatial  units  are  millimeters.  Low  thermal  conductance  (left) 
causes  the  heated  zone  to  spread  more  in  the  axial  direction,  making  the 
surface  temperature  distribution  broader  and  the  overall  temperature  lower. 
High  thermal  conductance  (right)  causes  the  heat  to  transfer  into  the  matrix 
near  the  irradiation  surface  so  that  the  overall  temperature  is  higher,  resulting 
in  a  sharper  gradient. 
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space,  so  no  temperature  maps  of  these  heating  geometries  are  complete. 


4.2  Scanning  Thermal  Microscopy 

In  order  to  verify  some  of  the  trends  in  interface  thermal  conductance  pre¬ 
dicted  by  the  MD  simulations,  experimental  approaches  were  developed  to 
measure  heat  flow  across  an  interface  directly.  Scanning  thermal  microscopy 
(SThM)  is  a  class  of  nanoscale  surface  thermal  measurements  that  employ 
a  surface  probe  with  some  thermal  functionality.  The  Anasys  Instruments 
NanoIR2  is  equipped  with  a  module  that  enables  the  relative  heat  flux  be¬ 
tween  a  special  tip,  shown  in  Figure  11,  and  the  sample. 


Figure  11:  SEM  of  Thermalever  probe  from  Anasys  Instruments 

(http://www.anasysinstruments.com/products/thermal-probes-tips-afm- 

sthm/) 

Typically,  SThM  is  used  to  generate  images  using  thermal  conductivity 
as  a  contrast  mechanism.  Attempts  to  use  probes  such  as  the  Thermalever 
probe  shown  in  Figure  11  are  complicated  by  the  many  conductive  pathways 
between  the  measurement  device  and  the  sample.  Careful  calibration  can 
be  performed,  but  it  requires  a  specialized  and  dedicated  apparatus,  usually 
under  high  vacuum  and  controlled  environment. 

In  practice,  the  SThM  probe  is  heated  with  a  bias  voltage  to  a  few  de¬ 
grees  above  ambient  temperature  (or  above  the  sample  temperature),  and 
the  voltage  across  the  probe  is  brought  to  zero  by  changing  the  value  of  a 
variable  resistor  placed  opposite  the  probe  in  a  wheatstone  bridge  configu¬ 
ration.  As  the  temperature  of  the  probe  changes,  so  does  its  resistance,  and 
the  SThM  signal  consists  of  the  voltage  that  appears.  Positive  voltages  result 
when  the  resistance  of  the  probe  increases  as  a  result  of  heat  not  flowing  out 


17 

Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 


of  the  probe,  and  negative  voltages  indicate  excess  heat  flowing  out  of  the 
probe.  Thus,  in  most  STliM  images  bright  regions  correspond  to  low  thermal 
conductance  and  dark  regions  correspond  to  high  thermal  conductance. 

Even  without  quantitative  calibration,  useful  information  can  still  be  ob¬ 
tained  from  the  probe,  but  only  relative  properties  between  multiple  samples 
can  be  reliably  determined.  By  subjecting  the  surface  to  a  constant  heat 
flux  along  a  prescribed  direction,  relative  changes  in  temperature  are  asso¬ 
ciated  with  changes  in  the  SThM  voltage.  These  changes  will  be  different 
for  different  materials,  so  the  temperature  gradient  across  two  dissimilar  ma¬ 
terials  is  unquantifiable  without  additional  knowledge  of  the  heat  transfer 
mechanisms  in  the  immediate  neighborhood.  Nevertheless,  for  two  samples 
of  similar  composition  but  different  interfaces  the  relative  interface  thermal 
conductance  can  be  calculated. 

When  the  sample  temperature  is  close  to  the  ambient  temperature,  so  that 
losses  to  the  environment  are  minimized,  it  can  be  shown  that  the  thermal 
resistance  R  across  an  interface  of  two  materials  having  different  effective 
thermal  conductivities  can  be  written  as 

R  =  m(Tl)(Sl-S2)  (2) 

where  Sn  is  the  SThM  signal  (voltage)  in  region  n,  and  m  is  a  smoothly 
increasing  function  of  temperature.  For  small  changes  in  temperature  (ie, 
small  transverse  heat  fluxes  as  shown  in  Figure  12)  m  is  approximately  linear. 
Thus,  the  relative  interface  resistance  between  two  interfaces  can  be  measured 
by  the  ratio  of  the  slopes  of  A S(Qt),  where  Qt  is  the  heat  flux  in  the  plane 
of  the  surface. 


Figure  12:  Schematic  for  semi-quantitative  SThM  measurements. 

Several  specimens  consisting  of  a  single  carbon  fiber  in  a  BMI  matrix 
were  fabricated,  and  the  samples  were  polished  to  expose  the  cross-section 
of  the  fiber  at  the  surface.  The  processing  of  test  specimens  for  the  SThM 
measurements  deviates  from  standard  procedures  for  neat  resin  and  prepreg 
composites.  Individual  carbon  fibers  are  suspended  or  draped  across  a  wide, 
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flat  silicone  mold.  Caution  must  be  taken  while  pouring  and  curing  the  resin 
to  avoid  trapped  bubbles  and  to  preserve  the  relative  position  and  orientation 
of  sparsely  distributed  fibers  through  buckling  of  the  neat  resin  during  cure. 
The  two  components  (one  a  liquid  resin,  the  other  a  dry  crystalline  powder) 
are  stirred  by  hand  in  a  hot  glycerin  bath  at  130°  C  until  a  uniform  melt 
is  obtained.  When  the  viscosity  decreases  during  melting  the  sheer  forces 
become  insufficient  to  break  up  the  dry  component,  and  the  mixture  must 
be  cooled  to  increase  the  viscosity.  This  is  repeated  several  times  until  a 
uniform  mixture  is  obtained. 

Once  the  components  are  mixed  the  vessel  is  placed  in  a  vacuum  oven 
preheated  to  150°  C  and  the  pressure  is  reduced  to  degas  the  mixture.  The 
resin  is  held  at  elevated  temperature  and  low  pressure  for  5  minutes  and  then 
removed  from  the  oven.  It  is  then  poured  carefully  into  the  mold  containing 
the  carbon  fibers  and  the  mold  is  placed  back  in  the  vacuum  oven  at  150°  C 
and  low  vacuum  for  30  min  more.  By  this  time  the  resin  has  started  to  gel. 
The  pressure  is  then  increased  to  atmospheric  pressure  and  the  temperature  is 
raised  to  the  cure  temperature  of  240°  C  for  4  hours.  In  order  to  minimize  the 
effect  of  topography,  samples  are  polished  according  to  the  recipe  in  Table  1 
until  the  difference  between  the  fiber  surface  and  the  surrounding  matrix  is 
less  than  200  nm.  Care  must  be  taken  not  to  polish  too  vigorously  or  too 
slowly  so  as  not  to  create  a  large  difference  between  the  fiber  and  matrix 
surfaces  because  of  the  difference  in  hardness  between  the  two  materials. 

A  variety  of  surface  preparation  techniques  were  attempted,  including  ion 
etching,  focused  ion  beam  milling,  chemical  etching,  and  microtome  cutting. 
However,  no  technique  consistently  produces  useable  samples  more  often  than 
the  traditional  polishing  technique  described  in  Table  1. 

Figure  13  shows  an  atomic  force  micrograph  (AFM)  of  one  such  sample. 
To  measure  Si  and  S2  precisely  requires  a  perfect  interface  boundary.  In 
practice,  however,  this  is  never  achieved.  Key  challenges  in  sample  prepara¬ 
tion  include  maintaining  a  minimum  height  difference  between  the  fiber  and 
the  matrix,  and  a  minimum  gap  between  the  fiber  and  the  matrix  at  the 
surface.  The  height  difference  between  the  fiber  and  the  sample  will  not  con¬ 
tribute  significantly  to  errors  in  the  STHM  signal  as  long  as  the  difference  is 
not  greater  than  one  third  the  height  of  the  probe.  To  account  for  a  physical 
gap  between  the  fiber  and  the  matrix  we  extrapolate  the  linear  signal  S(x ) 
where  x  is  the  position  along  the  trace.  This  approach  introduces  negligible 
error  as  long  as  the  gap  size  in  either  direction  is  less  than  the  probe  height. 
The  sample  shown  in  Figure  13  has  a  height  difference  of  less  than  200  nm 
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Table  1:  Polishing  procedure  for  the  preparation  of  samples  for  SThM  mea¬ 
surements. 

1.  Using  a  diamond  saw,  cut  slices  of  the  material  at  right  angles  to  the 
fiber  into  3mm  thick  slabs. 

2.  Attach  slabs  to  a  fixed  substrate  using  glue  or  crystal  bond. 

3.  Grind  the  top  surface  using  400  grit  SiC  paper  at  100  rpm  until  the 
surface  is  parallel  to  the  mount  and  uniformly  roughened.  Rinse  well 
with  DI  water. 

4.  Polish  at  successively  hirer  grits  (600,  800,  &  1200  grit),  pausing  every 
30-45  seconds  during  each  step  to  insure  that  the  surface  is  being  uni¬ 
formly  abraded.  If  large  streaks  develop  and  are  not  polished  easily,  go 
back  to  the  previous  step.  Clean  thoroughly  with  DI  water  after  each 
step. 

5.  Polish  with  0.3  /irn  alumina  suspension  on  ChemoMet  lapping  cloth  for 
5-10  minutes  to  remove  scratches. 

6.  Wash  with  detergent  (Microsoap  or  equivalent)  in  water,  then  rinse 
with  isopropanol  and  methanol/ethanol  to  remove  any  residue.  Dry 
under  flowing  nitrogen. 
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and  a  gap  width  of  less  than  500  nm,  which  is  sufficiently  small  for  the  probes 
that  were  used. 


Figure  13:  Topographical  map  of  a  cross-section  of  a  single  carbon  fiber 
emerging  from  the  surrounding  BMI  matrix  obtained  via  contact  mode  AFM. 

After  the  surface  is  prepared,  the  sample  is  mounted  on  a  nickel  puck 
for  placement  in  the  SThM  apparatus.  A  local  heat  source  along  one  edge 
is  used  to  generate  a  net  flux  along  the  plane  of  the  sample  surface,  and 
the  resulting  SThM  signal  is  analyzed  according  to  Equation  2  to  find  the 
relative  interface  thermal  resistance. 

Figure  14  shows  a  thermal  micrograph  for  one  sample  and  the  plots  of 
AS(Qt)  for  a  carbon  fiber  with  and  without  an  alumina  coating  as  a  thermal 
barrier.  The  obvious  increase  in  thermal  resistance  for  an  insulating  coating 
demonstrates  the  type  of  information  that  can  be  obtained  from  SThM  in 
this  way.  For  the  samples  in  Figure  14,  PAN-based  carbon  fibers  (Hexcel 
T300)  were  coated  with  200  nm  of  alumina  using  atomic  layer  deposition 
(ALD).  After  coating,  the  interfacial  thermal  resistance  between  the  fiber 
and  the  matrix  increased  by  a  factor  of  2.18.  This  change  is  consistent  with 
MD  predictions  for  the  effect  of  a  roughened  or  functionalized  surface,  so 
it  is  still  unclear  what  the  underlying  mechanism  is  for  the  increase.  At 
the  time  of  the  preparation  of  this  report,  investigations  are  underway  for 
samples  that  represent  the  types  of  surface  modifications  modeled  by  the 
MD  systems,  including  varying  the  surface  structure,  functionalization,  and 
presence  of  carbon  nanotubes. 
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Interfacial  SThM  Signal  Change 
With  Heat  Flow 


Figure  14:  Scanning  thermal  microscopy  of  carbon  fibers  in  a  BMI  matrix 
with  a  transverse  heat  flow,  (left)  SThM  Micrograph.  Inset  shows  the  varia¬ 
tion  of  the  signal  with  position  along  the  center  line,  (right)  Signal  differences 
at  the  interface  as  a  function  of  transverse  heat  flow.  The  higher  slope  for  the 
alumina-coated  sample  indicates  a  higher  thermal  resistance  at  the  interface. 


5  Conclusion 

An  ICMSE  approach  to  understanding  laser  interactions  in  OMCs  was  de¬ 
veloped  across  multiple  length  scales.  Molecular  dynamics  and  nanoscale 
surface  probes  were  used  together  to  determine  the  magnitude  and  influ¬ 
ence  of  interface  thermal  resistance  between  carbon  fibers  and  the  polymer 
host.  Finite  element  analysis  was  used  along  with  high  resolution  measures 
of  local  temperature  gradients  to  elucidate  the  influences  of  constituent  ma¬ 
terial  properties  on  optical  and  thermal  energy  transport  in  these  materials. 
As  part  of  this  effort,  a  novel  microscale  temperature  measurement  tech¬ 
nique  was  developed  using  gold  nanorods.  The  knowledge  and  techniques 
developed  during  this  project  will  be  invaluable  in  future  efforts  to  engineer 
adaptive  responses  to  laser  irradiation  at  the  fiber-matrix  interface. 
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Appendix  A.  Molecular  Dynamics  for  Modeling  of  In¬ 
terface  Thermal  Resistance 
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ABSTRACT:  The  rapid  heating  of  carbon-fiber-reinforced 
polymer  matrix  composites  leads  to  complex  thermophysical 
interactions  which  not  only  are  dependent  on  the  thermal 
properties  of  the  constituents  and  microstructure  but  are  also 
dependent  on  the  thermal  transport  between  the  fiber  and 
resin  interfaces.  Using  atomistic  molecular  dynamics  simu¬ 
lations,  the  thermal  conductance  across  the  interface  between  a 
carbon-fiber  near-surface  region  and  bismaleimide  monomer 
matrix  is  calculated  as  a  function  of  the  interface  and  bulk 
features  of  the  carbon  fiber.  The  surface  of  the  carbon  fiber  is 
modeled  as  sheets  of  graphitic  carbon  with  (a)  varying  degrees 
of  surface  functionality,  (b)  varying  defect  concentrations  in 
the  surface-carbon  model  (pure  graphitic  vs  partially  graphitic), 

(c)  varying  orientation  of  graphitic  carbon  at  the  interface,  (d)  varying  interface  saturation  (dangling  vs  saturated  bonds),  (e) 
varying  degrees  of  surface  roughness,  and  (f)  incorporating  high  conductive  fillers  (carbon  nanotubes)  at  the  interface.  After 
combining  separately  equilibrated  matrix  system  and  different  surface-carbon  models,  thermal  energy  exchange  is  investigated  in 
terms  of  interface  thermal  conductance  across  the  carbon  fiber  and  the  matrix.  It  is  observed  that  modifications  in  the  studied 
parameters  (a— f)  often  lead  to  significant  modulation  of  thermal  conductance  across  the  interface  and,  thus,  showcases  the  role 
of  interface  tailoring  and  surface-carbon  morphology  toward  thermal  energy  exchange.  More  importantly,  the  results  provide  key 
bounds  and  a  realistic  degree  of  variation  to  the  interface  thermal  conductance  values  at  fiber/matrix  interfaces  as  a  function  of 
different  surface-carbon  features. 

KEYWORDS:  carbon  fibers,  BMI  resin,  molecular  dynamics,  interfaces,  thermal  conductance 
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1.  INTRODUCTION 

Today,  laser  technology  is  being  used  in  a  wide  variety  of  fields 
including  communications,  industrial  and  environmental 
applications,  medicine,  and  materials’  research  and  develop¬ 
ment  (R&D).  Some  of  the  key  areas  in  materials’  R&D  include 
mass  spectral  analysis  using  matrix-assisted  laser  desorption/ 
ionization  (MALDI)1  and  fabrication  and  machining  of  micro-/ 
nano-patterned  surfaces  and  thin  films  using  pulsed  or 
continuous  lasers  (such  as  fabrication  of  nanoelectronic  devices 
and  patterning  of  polymeric  nanocomposites).2’3  Our  specific 
interest  lies  in  how  different  lasers  interact  with  the  polymeric 
composite  materials,  what  parameters  govern  the  interaction 
behavior,  and  how  these  laser— composite  interactions  lead  to 
evolution  of  composites’  thermophysical  properties  over  time. 

The  basic  mechanism  of  laser  interaction  with  materials 
involves  many  nonequilibrium  processes  (thermal,  mechanical, 
chemical,  and  physical)  at  varying  length  and  time  scales  caused 
by  fast  deposition  of  laser  energy  and  its  dissipation  in  various 
forms.4’'  Because  of  the  rapid  heating  rate,  these  non- 
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equilibrium  processes  cause  spatial  variability  in  thermophysical 
and  mechanical  properties  of  the  matrix  such  as  local  heating, 
local  matrix  degradation  due  to  sudden  local  increase  in 
temperature,  pressure  wave  formation,  etc.4’'  In  the  context  of 
polymeric  composites,  it  is  often  desired  to  delay  the  onset  of 
ablation  and  keep  the  composite  intact  by  tailoring  the 
interfaces  between  matrix  and  fillers  (such  as  carbon  fibers, 
carbon  nanotubes,  exfoliated  graphite/graphene,  and  metallic 
nano-/micro-particles,  etc.)  by  modulating  the  interactions 
between  the  two  components.  Since  the  ablation  onset  is 
primarily  driven  by  a  rise  in  temperature,  it  becomes  crucial  to 
understand  individual  thermal  conduction  characteristics  of  the 
matrix  and  fillers  (thermal  conductivity),  as  well  as  of  their 
mutual  interface  (interface  thermal  conductance). 
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Within  the  framework  of  carbon-fiber-based  nano-/micro- 
composites,  several  experimental  techniques  (laser  flash, 
photoacoustic,  3 ft),  thermoreflectance,  infrared  microscopy, 
etc.)  have  been  employed  to  measure  the  thermal  conductivity 
of  individual  components  as  well  as  the  "effective”  thermal 
conductivity  of  the  nanocomposites.6’  However,  because  of  the 
small  length  scales  involved,  it  has  been  significantly  challenging 
to  directly  and  accurately  measure  the  nanocomposite  interface 
features,  both  morphological  (roughness,  interfacial  bonding 
nature,  fractional  contact  area,  etc.)  and  thermal  transport 
characteristics.  Most  of  the  experimentally  reported  values  of 
interface  thermal  conductance  for  carbon-based  nanocompo¬ 
sites  are  evaluated  indirectly  (fitting  the  measured  experimental 
data  using  different  mathematical  models,  often  requiring 
thermal  properties  of  individual  components  as  input 
parameters).8-14 

On  the  contrary,  because  of  the  sharpness  of  molecular 
interfaces,  atomistic  molecular  dynamics  (MD)  simulations1' 
offer  an  alternative  as  well  as  a  direct  approach  to  calculate  the 
thermal  energy  exchange  across  matrix/filler  interfaces  in  terms 
of  interface  thermal  conductance.  In  these  simulations,  only  the 
details  of  microscopic  interactions  need  to  be  specified  (often 
extracted  using  quantum  chemical  calculations)  and  no  further 
assumptions  need  to  be  made  in  order  to  investigate  the 
equilibrium  or  nonequilibrium  phenomenon  happening  at 
these  time  scales.  To  date,  within  the  framework  of  carbon/ 
polymeric  interfaces,  most  molecular  modeling  reports  have 
estimated  thermal  conductance  primarily  with  either  carbon 
nanotubes  or  single/few-layer  graphene  (1—2  nm  fill¬ 
ers).14'16-'0  Furthermore,  in  many  of  these  studies,  interface 
functionalization  has  been  used  as  the  key  parameter  to  tailor 
the  interfacial  thermal  transport  characteristics.18’20'23’24'27'29’30 

In  principle,  apart  from  interface  functionalization,  several 
other  parameters  could  be  altered  to  modulate  the  thermal 
transport  across  carbon-fiber/matrix  interfaces.  This  could 
include  interface  characteristics  such  as  interface  roughness,  the 
nature  of  dangling  bonds,  incorporation  of  high  conductive 
fillers  (nanotubes,  metal  nanoparticles)  at  the  interface;  carbon- 
fiber  characteristics  such  as  its  amorphous  or  graphitic  nature 
and  graphitic  orientation  near  the  interface;  and  matrix  features 
such  as  degree  of  curing,  local  crystallinity,  or  molecular  order 
near  the  interface,  etc.  To  the  best  of  the  authors’  knowledge, 
these  parameters  have  either  not  been  explored  or  been 
explored  in  a  very  limited  way  from  the  perspective  of  interface 
thermal  conductance  across  carbon-fiber/matrix  interfaces. 

This  work  addresses  several  of  the  highlighted  issues  to  (a) 
model  different  types  of  carbon-fiber  surface  models  ranging 
between  ~5  and  ~17  nm  (mimicking  the  near-surface  region  of 
carbon  fiber  to  a  much  larger  scale  than  what  is  reported  to 
date);  (b)  model  high-temperature  BMI  monomeric  resins 
(specifically  Matrimid  BMI-5292  because  of  its  high  temper¬ 
ature  stability)  with  atomistic  resolution;  and  (c)  model  the 
effect  of  the  aforementioned  parameters  (including  functional¬ 
ization)  on  carbon-fiber/BMI-matrix  interfacial  thermal  trans¬ 
port  in  a  single  comprehensive  study. 

It  is  important  to  point  out  that  while  modeling  thermal 
transport  across  carbon-fiber/ crosslinked  resin  is  ideal,  the 
following  study  investigates  the  transport  characteristics  across 
its  uncrosslinked  counterpart  because  of  the  following  reasons. 
(A)  Cross-linking  reactions  in  BMI  are  notably  more  complex 
than  well-studied  epoxy-based  resins.  For  example,  unlike 
curing  in  epoxy,  BMI  curing  consists  of  five  different  reaction 
pathways;31'32  (B)  Specific  reaction  pathways  are  temperature 


dependent,  and  their  relative  importance  to  experimental 
conditions  are  not  well-quantified.  (C)  Inclusion  of  cross- 
linking  would  significantly  broaden  the  scope  of  the 
investigation,  i.e.,  investigating  the  effect  of  the  degree  of 
curing  and  the  effect  of  different  reactions  and  their  relative 
occurrence  on  investigated  properties,  etc.  We  do,  however, 
discuss  the  possible  modulation  of  interfacial  thermal  transport 
because  of  curing  phenomenon  toward  the  end  of  the  article. 
Understanding  thermal  transport  across  the  carbon-fiber  and 
monomeric  matrix  interface  is  also  pertinent  from  the 
perspective  of  laser-assisted  curing  of  carbon-fiber  composites. 
Here,  the  interfacial  thermal  transport  is  expected  to  govern  the 
heat  flow  from  carbon  fiber  to  the  matrix  (most  of  the  laser 
energy  is  absorbed  by  the  carbon  fiber),  subsequently 
determining  the  temperature  rise  and  distribution  in  the 
monomeric  matrix  and  thus  affecting  its  curing  kinetics.31'32 

2.  SIMULATION  SETUP  AND  METHODOLOGY 

The  section  is  divided  into  five  subsections,  dealing  with  the  (a) 
preparation  of  BMI  matrix,  (b)  preparation  of  different  surface-carbon 
and  (c)  interface  models,  (d)  preparation  of  combined  systems  for 
thermal  conductance  calculations,  and  (e)  simulation  methodology  for 
thermal  conductance  calculations. 

2.1.  Matrix  Preparation.  The  BMI  resin  (Matrimid  BMI-5292)  is 
composed  of  4,4'-bismaleimidodiphenyl  methane  (BMPM)  and  0,0'- 
diallylbisphenol  A  (DABPA)  monomers,  often  in  a  1:1  molar  ratio.33'34 
The  schematic  representations  of  these  monomers,  along  with  their 
molecular  structures,  are  shown  in  Figure  1.  As  a  part  of  the  initial 


(BMPM)  (DABPA) 


Figure  1.  Molecular  structure  and  schematic  representation  of  the  two 
monomer  components  of  the  Matrimid  BMI-5292  resin.  Color 
scheme:  cyan,  carbon;  blue,  nitrogen;  red,  oxygen;  silver,  hydrogen. 


simulation  setup,  these  monomers  were  constructed  in  Materials 
Studio.1^  Among  several  all  atom  force  fields  tested  (CVFF,  '6  PCFF,1 
and  GAFF38'39)  for  successful  generation  of  bonded  and  nonbonded 
parameters,  only  the  general  amber  force  field  (GAFF)  was  able  to 
generate  all  of  the  required  parameters.  Hence,  GAFF  parameters  were 
employed  to  model  the  BMI  matrix  as  well  as  all  carbon-surface 
models  for  this  study.  It  is  often  reported  that  when  employing  GAFF 
force  field,  restrained  electrostatic  potential  (RESP)  charges  are  best 
suited  to  model  organic  or  biomolecules.38’ 9  For  our  simulations, 
RESP  charge  estimations  and  fittings  were  performed  using 
GAUSSIAN09  software40  and  RESP  code  within  the  framework  of 
the  RE.D.  program.41  The  LAMMPS  package  was  employed  to  carry 
out  all  molecular  dynamics  simulations.42  For  van  der  Waals 
interactions,  a  distance  cutoff  value  of  12  A  was  used  while  PPPM 
(particle— particle— particle  mesh)  technique  was  used  to  model  long- 
range  Columbic  interactions  as  implemented  in  LAMMPS. 

First,  a  system  of  ~9500  atoms  was  created  by  repheating  both  BMI 
monomers  108  times  (total)  along  x-,  y-,  and  z-directions.  Following 
the  initial  minimization,  MD  simulations  were  carried  out  in  NVT 
(canonical)  ensemble  to  equilibrate  the  temperature  at  300  K  for  200 
ps,  followed  by  NPT  (isothermal— isobaric)  ensemble  to  equilibrate 
the  density  (and  pressure)  for  5  ns  with  independent  barostats  along 
the  x-,  y-,  and  z-directions.  Afterward,  200  ps  equilibration  was 
performed  in  NVE  (microcanonical)  ensemble  to  confirm  the  energy 
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(d)  Interface  Functionalization 


5  nm 


Wavelength 
2.5  nm 


1.7  nm 


Capped  Nanotube 


(e)  Interface  Roughness 


(g)  Interface  Saturation 


Figure  2.  Several  schematics  of  generated  carbon-fiber  surface  models  pertaining  to  (a)  graphitic  orientation;  (b)  incorporation  of  sp3  defects  within 
the  graphitic  structure  (structure  before  and  after  bonding  is  also  shown;  purple,  sp3  carbon  atoms;  cyan,  sp2  carbon  atoms);  (c)  effective  fiber 
thickness;  (d)  interface  functionalization  (close-up  interface  region  of  two  functionalized  cases  are  also  shown);  (e)  interface  roughness  with 
different  roughness  wavelengths  and  amplitude  (perspective  view  is  shown  to  reveal  two-dimensional  roughness  features);  (f)  incorporation  of 
bonded  CNTs  at  the  interface  (close-up  picture  of  nanotube  cap  is  also  shown);  and  (g)  saturation  of  interface  bonds  (saturated  H  atoms  are  shown 
as  blue  spheres  in  close-up  perspective). 
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conservation  of  the  equilibrated  system.  For  these  simulations, 
periodic  boundary  conditions  were  employed  in  all  three  orthogonal 
directions  to  mimic  bulk-like  behavior.  A  time  step  of  1.0  fs  was  used 
during  the  equilibration  phase.  After  successful  equilibration,  the 
system  was  subjected  to  a  series  of  additional  simulations  for 
estimation  and  experimental  comparison  of  the  glass  transition 
temperature,  room  temperature  density,  coefficient  of  thermal 
expansion,  and  thermal  conductivity  for  the  validation  of  the  employed 
force  field,  which  are  reported  elsewhere.43 

2.2.  Carbon-Surface  Models:  Bulk  Features  Preparation. 
Intrinsic  carbon-fiber  parameters  such  as  graphitic  order  and  its 
orientation  with  respect  to  heat  flow  direction,  content  of  amorphous 
carbon  within  the  carbon  fiber,  and  sp2/sp3  bond  ratio  are  important 
parameters  which  impact  its  bulk  thermal  transport  characteristics.  In 
principle,  these  parameters  could  also  affect  how  interfaces  of  such 
fibers  with  organic  matrices  respond  to  an  applied  thermal  gradient.  In 
order  to  explore  this,  we  generated  several  of  such  models 
incorporating  differences  in  the  aforementioned  bulk  carbon-fiber 
features  which  are  discussed  as  follows. 

2.2.7.  Graphitic  Orientation.  In  order  to  investigate  the  impact  of 
graphitic  layer  orientation  on  interfacial  thermal  transport,  several 
carbon-surface  models  were  created  with  several  different  orientation 
angles  of  graphitic  layers,  ranging  from  0°  to  90°  with  respect  to  the  re¬ 
direction  (parallel  to  the  interface).  We  will  refer  to  these  cases  as  Gq, 
where  6  represents  the  orientation  angle  of  the  graphitic  layers  with 
respect  to  the  interface.  Several  models,  with  orientation  angles  of  0°, 
8°,  16°,  24.5°,  33°,  41°,  51°,  66°,  and  90°  were  created  using  in-house 
scripts  with  a  fiber  surface  thickness  of  10  nm.  The  details  of  building 
the  carbon-surface  models  with  different  orientations  are  outlined  in 
Supporting  Information  section  SI.  Figure  2a  shows  a  few 
representative  cases  of  carbon-surface  models  with  graphitic 
orientations  of  Gq)  c33,  and  C90. 

2.2.2.  Incorporation  of  Defects.  From  the  perspective  of  defects 
inclusion,  pristine  graphitic  and  amorphous  carbon  fibers  represent 
two  extremes  of  structural  order.  While  models  of  pristine  graphitic 
carbon  are  discussed  in  earlier  text,  interface  thermal  conductance 
across  amorphous  carbon  and  BMI  matrix  was  recently  modeled  by  us 
and  is  discussed  elsewhere.4 1  In  this  study,  we  create  partially  defective 
graphitic  carbon-fiber  surface  models  and  investigate  the  effect  of 
structural  graphitic  defects  within  carbon  fibers  on  thermal 
conductance  across  fiber/matrix  interfaces. 

The  defective  graphitic  carbon  models  can  either  be  generated  by 
incorporating  defects  in  the  pristine  graphite  or  incorporating  graphitic 
order  starting  from  fully  amorphous  carbon.44  For  this  study,  the 
former  route  is  chosen  to  create  graphitic  models  with  defects.4^  In 
this  approach,  we  create  interlayer  bonding  between  carbon  atoms  and 
modify  their  hybridization  from  sp2  to  sp3.  Here,  we  make  use  of  the 
fact  that,  in  AB  layer  stacking  of  graphite,  there  is  one  set  of  atom  pair 
which  lie  on  top  of  each  other  (looking  along  the  c-axis),  where  atoms 
are  3.4  A  apart.46  We  randomly  choose  a  predetermined  number  of 
such  pairs  (based  on  the  final  defect  concentration),  move  atoms  in 
each  pair  closer  toward  each  other  (l  A  each),  and  create  a  bond 
between  the  two.  The  bond  creation  is  followed  by  updating  the 
surrounding  topology  information  (adding  several  new  bond  angles 
and  dihedral  angles  and  deleting  two  improper  angles).  By  doing  so, 
several  models  were  created  with  volumetric  defect  concentrations 
ranging  up  to  ~6  defects/ nm3  for  three  graphitic  orientations,  namely, 
C0,  C41,  and  C90,  and  are  shown  in  Figure  2b,  which  also  shows  the 
structural  changes  before  and  after  the  bond  creation  for  a 
representative  case. 

2.2.3.  Fiber  Surface  Morphology  Thickness.  In  order  to  explore  the 
effect  of  thickness  of  fiber  surface  morphology  (effective  system  size  in 
MD  simulations),  several  models  were  created  with  varying  carbon- 
fiber  thicknesses  along  the  z-direction.  Here,  the  Qi  graphitic 
orientation  model  with  thickness  values  ranging  from  6  to  18  nm 
were  created  by  deleting  atoms  on  the  boundaries  of  a  much  longer 
C4i  model  along  the  z-direction.  Two  of  the  created  systems  are  shown 
in  Figure  2c. 

2.3.  Carbon-Surface  Models:  Interface  Features  Preparation. 

In  addition  to  bulk  carbon-fiber  surface  characteristics,  interface 


modification  can  also  provide  substantial  modulation  of  interfacial 
thermal  transport  characteristics  as  evident  from  previous  studies, 
mainly  focused  on  functionalization.18’20'23'24'27'29,30  Below,  we  detail 
other  interface  modifications  (in  addition  to  functionalization)  that  can 
be  employed  to  modify  interface  thermal  conductance  across  fiber/ 
matrix  interfaces. 

2.3.7.  Functionalization.  In  order  to  investigate  the  effect  of 
interface  functionalization,  five  different  carbon-surface  models  were 
created  using  graphitic  orientation  model  C90  with  a  surface 
functionalization  density  ranging  from  0.05  to  0.55  molecules/nm2. 
Among  different  graphitic  orientation  models,  the  Q>o  model  was 
chosen  for  interface  functionalization  studies  because  of  two  primary 
reasons;  (a)  in  literature,  most  investigations  dealing  with  the  effect  of 
interface  functionalization  on  thermal  conductance  have  been 
performed  on  C0  type  interfaces,  such  as  the  side  wall  of  CNTs  or 
surface  modification  in  single-layer/few-layer  graphene,  and  their 
outcomes  are  well-investigated;  and  (b)  the  Ggo  model  provides  the 
most  number  of  possible  functionalization  sites  because  of  the  most 
number  of  edge  carbon  atoms  among  nonzero  graphitic  orientation 
models.  For  the  C90  model,  interface  functionalization  was  carried  out 
by  bonding  one  of  the  monomers  (DAB PA)  to  the  edge  carbon  atoms 
of  a  C90  carbon  surface.  In  order  to  carry  out  bonding  at  the  interface, 
first,  a  combined  system  of  a  BMI-matrix  and  G90  graphitic  model  was 
created.  Then,  a  predetermined  number  of  DABPA  molecules  were 
slowly  dragged  to  the  carbon-fiber  surface  over  the  course  of  an  NVT 
MD  simulation.  Thereafter,  bonds  were  created  between  the  “hydroxyl 
O”  (of  DABPA)  and  “nearest  edge  carbon  atom”  using  an  in-house 
script,  also  updating  new  topology  information  (angles,  dihedrals,  and 
new  force-field  parameter  assignment,  etc.).  The  full  procedure  for 
functionalization  and  updating  topology  is  outlined  in  detail  in 
Supporting  Information  section  S2.  A  few  representative  function¬ 
alized  DABPA  monomers  (surrounding  matrix  removed  for 
clarification)  are  shown  in  Figure  2d. 

2.3.2.  Roughness.  The  effect  of  carbon-fiber  surface  roughness  on 
interfacial  thermal  transport  was  modeled  by  creating  rough  graphitic 
edges  for  a  Qo  graphitic  system  (Please  see  Supporting  Information 
section  S3  for  specific  details  regarding  generation  of  various 
roughness  models).  Here,  nine  carbon-surface  models  were  created 
with  varying  sinusoidal  roughness  features  (amplitude  and  wavelength) 
along  both  x-  and  y-directions.  Specifically,  roughness  with  three 
different  amplitudes,  i.e.,  1,  2  and  3  nm,  and  3  different  wavelengths, 
i.e.,  ~1.7,  ~2.5,  and  ~5  nm,  were  modeled.  A  few  schematic  carbon- 
surface  models  with  varying  roughness  features  are  shown  in  Figure  2e. 

2.3.3.  CNT  Incorporation  at  the  Interface.  In  addition  to 
functionalization  and  roughness,  incorporating  highly  conductive 
nanofillers  at  the  interface  can  also  significantly  modulate  interfacial 
thermal  transport.  For  example,  CNTs  grown  directly  in  carbon  fibers 
(also  known  as  fuzzy  carbon  fibers)  could  provide  additional  surface 
area  to  interact  with  the  matrix,  potentially  enhancing  overall  thermal 
transport  characteristics.  In  order  to  mimic  such  an  interface,  several 
models  were  created  by  covalently  bonding  ( 10,10)  end-capped  single¬ 
wall  CNTs  of  different  lengths  (~3.8,  ~6,  and  ~9  nm)  to  a  G90 
carbon-surface  model.  The  end-capping  of  the  CNTs  was  achieved 
using  NanoCap  software.47  The  details  of  the  bonding,  such  as  how 
the  junction  was  created  along  with  other  specific  details,  are  outlined 
in  Supporting  Information  section  S4.  The  schematic  representations 
of  all  three  CNT  bonded  carbon-surface  models  are  shown  in  Figure 
2f. 

2.3.4.  Interface  Saturation.  From  previous  discussion,  it  should  be 
noted  that  our  graphitic  carbon  surface  models  (C8—  C90)  have 
unsaturated  or  undersaturated  valence  for  the  edge  atoms.  In  order  to 
estimate  if  saturation  of  such  valence  could  have  significant  impact  on 
interfacial  thermal  transport,  we  saturated  four  of  the  nine  graphitic 
models,  namely,  C16,  C33,  C51,  and  C90,  with  hydrogen  atoms  for 
interface  thermal  conductance  calculations.  The  visual  differences 
between  unsaturated  and  saturated  surface-carbon  atoms  are  shown  in 
Figure  2g  for  the  C33  graphitic  model. 

Afterward,  each  of  the  prepared  carbon-surface  models  was 
subjected  to  equilibration  in  NVT,  NPT,  and  NVE  ensembles  for 
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200,  2000,  and  200  ps  with  time  steps  of  1,  1,  and  0.5  fs,  respectively, 
prior  to  combining  with  the  matrix. 

2.4.  Combined  Preparation.  To  calculate  interface  thermal 
conductance  across  the  fiber/ matrix  interface,  combined  systems  were 
created  using  relaxed  BMI-matrix  and  separately  equilibrated  surface- 
carbon  models  (graphitic  orientation,  defects,  fiber  thickness,  rough¬ 
ness,  and  interface  saturation)  to  mimic  the  fiber/matrix  interface 
(Figure  3a).  We  should  point  out  that  functionalized  and  CNT- 


Figure  3.  (a)  Schematic  representation  of  how  the  combined  systems 
were  prepared.  It  should  be  noted  that  the  equilibrated  structure 
(extreme  left)  has  a  lesser  volume,  (b)  Schematic  of  nonequilibrium 
molecular  dynamics  simulations  as  employed  in  this  investigation  for  a 
representative  case.  Color  scheme:  black,  fixed  atoms;  red,  hot 
thermostated  atoms;  blue,  cold  thermostated  atoms.  Thermal  energy  is 
directed  from  the  graphitic  carbon  surface  toward  the  matrix. 


incorporated  systems  were  combined  during  carbon-surface  model 
preparations  as  discussed  in  Supporting  Information  sections  S2  and 
S4,  respectively.  For  the  remaining  cases,  an  in-house  script  was 
written  to  merge  the  two  structures  in  order  to  generate  the  initial 
combined  geometries.  Within  this  code,  the  BMI  matrix  in  the 
equilibrated  slab-like  geometry  (4  times  replication  along  the  z- 
direction  of  the  original  system  followed  by  its  equilibration)  was 
combined  with  the  equilibrated  carbon-surface  models  to  generate 
initial  coordinates  for  modeling  matrix/fiber  interface  systems  in 
LAMMPS  format  (Figure  3a)  (~60,000  atoms  in  each  simulation). 
While  combining  both  structures,  a  5  A  gap  was  purposely  included  to 
avoid  overlap  of  atoms.  Interaction  between  the  matrix  and  fiber  was 
evaluated  using  van  der  Waals  interactions.  Prior  to  thermal 
conductance  simulations,  all  combined  systems  were  equilibrated  for 
200,  2000,  and  200  ps  in  NVT,  NPT,  and  NVE  ensembles,  respectively, 
with  respective  time  steps  of  1,  1,  and  0.5  fs.  We  should  point  out  that, 
during  equilibration,  the  orientation  angle  for  C90  cases  was  eventually 
equilibrated  to  ~83°,  while  orientation  angles  for  other  graphitic  cases 
did  not  change  to  a  noticeable  degree. 

2.5.  Interface  Thermal  Conductance  Calculations.  The 
thermal  conductance  across  BMI-matrix  surface-carbon  models’ 
interfaces  was  calculated  using  nonequilibrium  molecular  dynamics 
(NEMD)  simulations  from  a  slab-like  geometry  based  on  Fourier  law 
approach48  and  is  also  schematically  depicted  in  Figure  3b.  As  seen 
from  the  figure,  the  outer  edges  of  the  elongated  slab  (l  nm)  were 
fixed,  while  the  adjacent  1.5  nm  bins  were  employed  as  hot  and  cold 
thermostats,  where  a  temperature  rescaling  algorithm  was  employed  to 
keep  the  respective  temperatures  to  350  and  250  K.  NVE  ensemble 
was  employed  for  nonthermostated  regions.  In  order  to  keep  the 
regions  at  their  specified  temperatures,  energy  was  continuously  added 
and  taken  off  from  the  hot  and  cold  regions,  respectively.  Our  previous 


studies  on  silicon  and  polymer  composites  have  shown  that  thermal 
conductivity  and  interface  thermal  conductance  prediction  is  not 
noticeably  affected  by  choice  of  thermostats  (Nose— Hoover, 
Berendsen,  and  temperature  rescaling).49-51 

By  applying  thermostats,  a  temperature  gradient  is  established  over 
time  along  the  elongated  slab  direction,  leading  to  a  steady-state 
temperature  profile.  Furthermore,  near  the  interface,  a  temperature 
discontinuity  develops  due  to  a  mismatch  between  the  vibrational 
properties  of  the  two  media.  The  temperature  profile  can  then  be 
calculated  by  dividing  the  slab  into  a  predefined  number  of  thin  slabs 
with  equal  thickness  (2  A  for  our  case)  and  calculating  the  temperature 
of  each  bin  using  the  following  equation. 


T;  = 


£ 


(i) 


where  N,  is  the  number  of  atoms  in  the  zth  slab  and  m*.  and  v 
correspond  to  atomic  mass  and  velocity  of  atom  k,  respectively.  For 
the  current  study,  each  of  the  2  A  bins  contained  ~500  atoms.  As  an 
example,  one  such  steady-state  temperature  profile  for  a  representative 
case  is  shown  in  Supporting  Information  section  S5  along  with  the 
calculation  of  temperature  discontinuity,  AT.  The  interface  thermal 
conductance  was  then  calculated  using  the  following  equation. 


A  = 


d 

AAT 


(2) 


where  A,  A,  and  Q  respectively  correspond  to  interface  thermal 
conductance,  cross-sectional  area,  and  heat  flux  as  estimated  by  the 
energy  dumped  into  hot  (taken  out  from  cold)  thermostats.'-  The 
heat  flux  calculation  is  further  discussed  in  Supporting  Information 
section  S6  using  the  cumulative  energy  dumped  (or  taken  out) 
estimates  for  a  representative  case.  For  all  cases,  initial  simulations 
were  run  for  about  1  ns  in  order  to  achieve  steady-state.  Following 
that,  simulations  were  further  run  for  4  ns  for  data  collection  for 
temperature  gradient  and  heat  flux  calculations. 


3.  RESULTS  AND  DISCUSSION 

Next,  we  present  and  discuss  the  effect  of  various  parameters  on 
thermal  transport  across  the  carbon-fiber/matrix  interface. 
Here,  unlike  the  previous  section,  we  have  not  divided  the 
discussion  between  carbon-surface  and  interface  models  but 
have  attempted  to  present  the  effect  of  different  build  models  in 
a  more  coherent  way.  We  wrap  up  this  section  with  a  few 
comments  on  possible  modifications  to  interface  thermal 
conductance  due  to  matrix  cross-linking,  which  is  not 
investigated  in  this  study. 

3.1 .  Effect  of  Functionalization.  We  start  this  section  with 
discussing  the  effect  of  surface  functionalization  density  on 
interface  thermal  conductance  as  shown  in  Figure  4.  The  figure 
clearly  shows  that  the  interface  functionalization  significantly 
increases  the  thermal  conductance  across  the  interface  with 
respect  to  unfunctionalized  interface  (C90),  increasing  up  to  4- 
fold  for  surface  functionalization  density  of  0.55  molecules/ 
nm2.  The  increasing  trend  and  relative  increase  is  in  excellent 
agreement  with  the  previous  literature.18’20’23'24’27’29'30  Here, 
the  increase  in  thermal  conductance  is  attributed  to  the 
interface  bonding  which  provides  additional  pathways  via 
covalent  bonds,  in  addition  to  van  der  Waals  interactions,  for 
better  thermal  energy  transport  across  the  interface. 

The  surface  functionalization  is  expected  to  occur  via  "end 
atoms”  of  polymeric  fragments,  whose  volumetric  concen¬ 
tration  for  oligomeric  or  polymeric  chains  is  expected  to  be  <SCl 
nm3.  In  this  context,  we  believe  that  the  largest  modeled  surface 
functionalization  density  of  0.55  molecules/nm2  is  noticeably 
higher  than  experimentally  achievable  surface  functionalization 
densities  of  attaching  amorphous  polymeric  fragments  on 
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Surface  Functionalization  Density  (nm‘2) 

Figure  4.  Plot  of  interface  thermal  conductance  as  a  function  of  surface 
functionalization  density  for  studied  functionalized  cases.  Two 
schematics  of  functionalized  interface  are  also  shown.  Relative  increase 
in  conductance  with  respect  to  unfunctionalized  C90  case  is  plotted  as 
well. 


carbon-fiber  surfaces.  Thus,  current  study  provides  an  estimate 
of  conductance  enhancement  bounds,  achievable  in  realistic 
experimental  conditions. 

3.2.  Effect  of  Interface  Roughness.  Generally,  when 
dealing  with  rough  interfaces,  there  is  an  overlap  region  where 
both  components  (carbon-fiber  surface  and  matrix  in  our  case) 
are  present.  Estimating  thermal  transport  across  such  cases 
introduces  the  concept  of  “interphase”  as  the  definition  of 
“interface”  or  a  distinct  “plane”  becomes  blurry.  Such  behavior 
is  depicted  in  Figure  5  for  one  of  our  model  systems  which 
shows  the  presence  of  both  graphitic  carbon  and  matrix  overlap 
over  a  3  nm  region.  In  principle,  one  could  define  various 
“planes”  or  “interfaces”  (as  depicted  in  Figure  5),  but  such  a 


Distance  along  slab  (nm) 

Figure  S.  Plot  of  steady-state  temperature  profile  across  rough 
interfaces  along  the  slab  direction.  Three  different  interfaces  are 
depicted  within  the  rough  interphase  with  dashed— dotted  lines  and  are 
discussed  in  Supporting  Information  section  S7.  “A”  and  “B”  denote 
opposite  boundaries  of  the  interphase.  Only  the  upper  halves  of  the 
error  bars  are  shown  for  clarity  purposes. 


definition  is  arbitrary  and  leads  to  different  values  of  interface 
thermal  conductance  based  on  the  chosen  interface  (please  see 
Supporting  Information  section  S7  for  further  discussion).  In 
order  to  avoid  that  ambiguity,  we  calculate  interphase  thermal 
conductance  between  points  A  and  B,  which  denote  the 
boundaries  of  the  interphase  region,  as  shown  in  Figure  5.  The 
resulting  interphase  thermal  conductance  is  plotted  in  Figure  6 
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Interphase  Thickness  (nm) 

Figure  6.  Plot  of  interphase  thermal  conductance  as  a  function  of 
interphase  thickness  for  three  different  roughness  wavelengths  of  5 
(red),  2.5  (green),  and  1.7  nm  (blue).  Black  data  points  show 
“effective”  thermal  conductance  across  the  same  width  for  the  flat  C90 
graphitic  model.  Only  the  upper  halves  of  the  error  bars  are  shown  for 
clarity  purposes. 

as  a  function  of  interphase  width  for  three  different  wavelength 
periodicities  and  compares  the  modulation  with  respect  to  the 
thermal  conductance  of  a  flat  interface  (C90)  for  similar  width 
(black  data  points). 

For  a  constant  interface  width,  the  data  show  that  the 
thermal  conductance  increases  with  decreasing  periodicity.  This 
increase  can  be  attributed  to  an  increase  in  carbon-fiber  surface 
area  which  interacts  with  the  matrix,  leading  to  enhanced 
thermal  conduction.  For  modeled  interface  widths,  our 
simulations  show  a  relative  increase  by  a  factor  of  2—3  in 
thermal  conductance  with  respect  to  their  flat  counterpart.  For 
a  constant  wavelength,  however,  the  data  show  that  the 
interphase  thermal  conductance  does  not  alter  significantly  and 
does  not  show  a  clear  trend  with  increase  in  interphase  width 
(or  roughness  amplitude).  While  increase  in  interphase  width 
increases  total  thermal  energy  exchange  (increase  in  surface 
area),  it  also  increases  the  width  over  which  the  thermal 
conductance  is  estimated.  These  two  features  compete  against 
each  other  in  interphase  thermal  conductance  estimation 
(increase  in  both  numerator  and  denominator),  leading  to  an 
overall  minor  variation  in  Figure  6.  Nevertheless,  it  is  important 
to  note  the  relative  increase  in  thermal  conductance  with 
respect  to  flat  interfaces  as  a  function  of  roughness  amplitude  as 
well  as  periodicity  of  roughness  undulation. 

3.3.  Effect  of  CNT  Incorporation  at  the  Interface. 
Similar  to  roughness,  incorporation  of  CNTs  at  the  carbon- 
fiber/matrix  interface  reignites  the  conductance  discussion  in 
terms  of  interphase.  In  this  context,  Figure  7  (inset)  shows  the 
estimated  thermal  conductance  across  the  interphase  as 
depicted  by  points  A  (fuzed  junction)  and  B  (nanotube  cap). 


26679  DOI:  10.1021/acsami.5b08591 

29  ACS  Appl.  Mater.  Interfaces  2015,  7,  26674-26683 


Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 


ACS  Applied  Materials  &  Interfaces 


Figure  7.  Plot  of  interphase  thermal  conductivity  as  a  function  of 
interphase  thickness  (or  CNT  length).  Blue  data  points  show  thermal 
conductivity  of  CNT-bonded  interphases  for  different  widths  while  red 
data  points  show  “effective”  thermal  conductivity  for  the  same  widths 
for  the  flat  C90  graphitic  model  with  no  CNTs.  Inset:  Plot  of 
interphase  thermal  conductance  as  a  function  of  interphase  thickness 
for  four  different  studied  cases:  flat  interface  (red);  3.8  nm  CNT 
(green);  6  nm  CNT  (cyan);  9  nm  CNT  (pink). 


Over  a  length  of  9  nm,  only  a  slight  decrease  in  conductance  is 
observed  which  is  indicative  of  a  more  conductive  interphase 
and  is  attributed  to  an  increase  in  effective  surface  area  for 
thermal  energy  transport."  '  The  ramifications  of  this  slight 
decrease  are  further  substantiated  in  Figure  7  (main)  which 
plots  the  interphase  thermal  conductivity  as  a  function  of 
interphase  thickness.  Also  shown  in  the  same  plot  is  the 
estimated  thermal  conductivity  of  the  same  region  (including 
the  interface)  if  only  the  matrix  were  present.  In  addition  to  the 
increase  in  effective  thermal  conductivity  of  the  interface  with 
CNT  length,  the  figure  shows  that  this  increase  is  also 
noticeable  (up  to  2-fold  for  modeled  cases)  with  respect  to  the 
effective  thermal  conductivity  of  the  pure  matrix  in  that  region. 
Based  on  the  observed  trends,  it  is  expected  that  higher  surface 
density  and  a  longer  length  of  covalently  bonded  CNTs  at 
carbon-fiber  surfaces  can  greatly  enhance  thermal  transport 
properties  across  the  fiber/matrix  interface. 

3.4.  Effect  of  Fiber  Surface  Morphology  Thickness.  For 
high  thermal  conductivity  materials,  such  as  CNTs  and 
graphite,  it  has  repeatedly  been  shown  that  the  system  size  of 
MD  simulations  affects  the  prediction  of  thermal  transport 
properties  because  of  their  extremely  large  phonon  mean  free 
paths.46'49'54’55  For  our  system,  modeled  graphitic  surfaces  of  10 
nm  also  present  a  similar  scenario.  In  order  to  investigate  if  the 
thickness  of  graphitic  surface  models  (equivalent  to  mimicking 
the  thickness  of  carbon-fiber  surface  morphology)  has  a 
consequential  effect  on  interfacial  thermal  transport,  we 
modeled  several  graphitic  lengths  for  the  C33  model  system. 
In  this  context,  Figure  8  shows  the  calculated  thermal 
conductance  values  for  different  studied  lengths  ranging  from 
~5  to  ~17  nm.  The  figure  clearly  shows  that  interface  thermal 
conductance  does  increase  significantly  with  increasing  graph¬ 
itic  length.  The  primary  reason  behind  such  an  increase  is 
attributed  to  the  inclusion  of  longer  wavelength  phonons  within 
the  graphitic  model  with  increasing  thickness,  which  can  couple 
with  the  BMI  matrix  when  they  propagate  to  the  interface. 
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Figure  8.  Plot  of  interface  thermal  conductance  as  a  function  of 
effective  fiber  thickness  for  the  C33  graphitic  model  case. 

While  the  simulations  were  performed  for  10  s  of  nm  of 
graphitic  thickness,  the  conductance  values  are  expected  to 
saturate  eventually,  when  fiber  dimensions  exceed  character¬ 
istics  mean  free  path  of  in-plane  phonons  in  graphite  (~^m). 

3.5.  Effect  of  Graphitic  Orientation.  Next,  we  discuss 
how  the  orientation  of  graphitic  layers  can  modulate  the 
thermal  energy  transport  across  the  interphase.  In  this  context, 
Figure  9  plots  and  shows  the  monotonic  decrease  of  the 
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Figure  9.  Plot  of  interface  thermal  conductance  as  a  function  of 
graphitic  orientation  for  unsaturated  (blue)  and  saturated  (red) 
interface  bonds.  The  numbers  next  to  the  red  data  points  show  the 
number  of  graphene  layer  edges  at  the  interface.  As  the  graphitic 
orientation  angle  increases,  so  do  the  number  of  graphene  layer  edges. 
Two  schematics  for  Q}3  and  Qo  graphitic  models  demonstrate  the  said 
increase  in  graphene  layer  edges. 


estimated  interface  thermal  conductance  as  a  function  of 
graphitic  orientation  angle,  9.  For  clearer  interpretation,  we 
separate  our  discussion  into  0°  (C0)  and  non-0°  (C8— C90) 
graphitic  orientation  cases.  The  reason  being  that  for  C8— C90 
cases,  graphene  layers  interacting  with  the  matrix  are  directly 
connected  to  the  hot  thermostat,  while  for  the  C0  case,  energy 
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has  to  propagate  from  the  hot  thermostated  graphene  layers  to 
the  end-graphene  surface  directly  interacting  with  the  matrix 
through  in-between  mediating  graphene  layers  via  van  der 
Waals  interactions  (See  Supporting  Information  section  S8  for 
better  pictorial  visualization). 

First,  the  exceptionally  large  value  (183  MWm-2  K-1)  of 
interface  thermal  conductance  for  case  C0  is  attributed  to  soft 
low-frequency  modes  propagating  across  graphene  layers. 
These  low-frequency  modes  can  couple  to  a  greater  extent 
with  the  low-frequency  vibrational  modes  of  the  matrix  leading 
to  enhanced  thermal  transport.  Recent  experiments  on  radial 
thermal  transport  across  large  multiwall  CNTs  have  suggested 
the  phonon  mean  free  path  across  the  out-of-plane  direction  of 
graphite  to  be  >100  nm.56  In  view  of  that,  it  is  expected  that  the 
predicted  conductance,  based  on  10  nm  graphitic  thickness,  is 
expected  to  increase  notably  for  much  thicker  flat  graphitic 
(~0°)  carbon-fiber  surfaces. 

On  the  other  extreme,  thermal  conductance  for  perpendic¬ 
ular  orientation  (C90)  was  observed  to  be  ~5-fold  lower  (~40 
MWm-2  K_1)  for  the  same  model  thickness  of  ~10  nm.  The 
low  value  is  attributed  to  the  small  thickness  (10  nm)  as  well  as 
relatively  stiff  (high-frequency)  propagating  modes.  When  such 
phonons  reach  the  interface,  the  probability  of  their  reflection 
at  the  interface  is  greater  than  soft  out-of-plane  modes  of 
graphite  (in  the  0°  case),  leading  to  a  lower  value  of  thermal 
conductance.  Figure  9  also  shows  that  interface  thermal 
conductance  is  increasing  with  decreasing  0  for  C8— C90  cases. 
This  increase  in  conductance  is  attributed  to  an  increase  in  the 
“effective  length”  of  the  graphite  (L/sin(0)).  Please  see 
Supporting  Information  section  S8  for  visualization  of  the 
“effective  length”.  In  these  cases,  phonons  predominantly 
propagate  along  graphitic  planes  (in-plane  modes)  due  to  the 
extremely  high  thermal  anisotropy  of  graphite  (>100  times) 
before  reaching  graphene  edges,  which  eventually  dictate  the 
thermal  transport. 

3.6.  Effect  of  Interface  Bond  Saturation.  Figure  9  also 
shows  the  interface  thermal  conductance  for  four  different 
graphitic  orientations  (C16,  C33,  C55,  and  C90)  where  we 
saturated  the  edge  carbon  atoms  by  hydrogen  atoms  (H 
atoms).  First,  for  constant  graphitic  orientation,  the  figure 
shows  that  saturated  models  always  show  higher  thermal 
conductance  than  their  unsaturated  counterpart.  The  added  H 
atoms  contribute  toward  enhancing  van  der  Waals  interactions, 
while  taking  minimal  volume  (small  a  and  short  C— H  bond 
length),  resulting  in  a  relative  increase  of  conductance. 

On  the  contrary,  unlike  unsaturated  cases,  we  see  an  increase 
in  conductance  with  increase  in  graphitic  orientation  angle.  We 
attribute  this  increasing  trend  to  the  increasing  number  density 
of  interacting  atoms  (H  atoms)  with  the  matrix  near  the 
interface.  As  noted  in  the  figure,  the  number  of  interacting 
graphene  layers’  edges  increase  with  an  increase  in  orientation 
angle,  which  in  turn  also  increases  the  surface  number  density 
of  H  atoms,  leading  to  enhanced  van  der  Waals  interactions 
(the  interactions  scale  as  N2,  N  being  the  number  of  interacting 
atoms).  Despite  the  shorter  effective  length  for  large  0  models, 
enhanced  interactions  dominate  and  lead  to  an  increase  in 
conductance,  also  increasing  the  respective  gap  between 
saturated  and  unsaturated  cases.  We  should  note  that  while 
this  is  a  somewhat  hypothetical  set  of  simulations,  it  does 
provide  important  information  on  how  the  number  density  of 
surface  interacting  atoms  can  modulate  interface  thermal 
conductance. 


3.7.  Effect  of  Graphitic  Defects.  Similar  to  graphitic 
orientation,  the  defects  within  carbon-fiber  surface  models  can 
also  affect  interface  thermal  conductance  as  shown  in  Figure  10. 
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Figure  10.  Plot  of  interface  thermal  conductance  as  a  function  of  sp' 
defect  concentrations  for  three  different  studied  graphitic  models,  C0, 
C41,  and  C90.  Schematics  of  these  models  for  ~3%  defect 
concentrations  are  shown  for  further  visual  clarity.  Furthermore, 
only  the  upper  halves  of  the  error  bars  are  shown  for  clarity  purposes. 
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The  figure  shows  significant  decrease  (~60%)  in  conductance 
for  the  C0  model  with  respect  to  modeled  defect  concen¬ 
trations,  while  the  decrease  is  relatively  minor  for  C41  (~20%) 
and  C90  (~0%)  graphitic  models.  For  the  former,  the  notable 
decrease  is  attributed  to  increased  perturbation  in  coherency  of 
out-of-plane  modes  with  increased  defect  concentrations  as 
these  defects  (sp3)  significantly  reduce  the  phonon  mean  free 
path  along  the  out-of-plane  direction.  However,  for  latter 
models  (C41  and  C90),  the  minor  or  negligible  decrease  is 
associated  with  the  observation  that,  in  these  cases,  sp3  bonds 
are  orientated  perpendicular  to  the  dominant  thermal  transport 
direction  (along  graphene  direction),  minimally  perturbing  the 
energy  carrying  propagating  phonons. 

3.8.  Effect  of  Matrix  Cross-Linking  on  Interface 
Thermal  Conductance.  Similar  to  intrinsic  fiber  properties 
and  interface  features,  it  is  expected  that  matrix  properties 
(cured  vs  uncured)  would  also  affect  interface  thermal 
conductance.  However,  this  parameter  space,  i.e.,  effect  of 
curing,  is  not  explored  in  this  study.  Unlike  modeling  of  an 
epoxy-based  matrix  curing  process,  which  often  requires 
modeling  a  single  amine— epoxy  reaction,5"  modeling  the 
curing  process  of  BMI-S292  resin  involves  five  separate 
reactions  (ene  reaction,  maleimide  homopolymerization, 
alternating  maleimide— adduct  co-polymerization,  ene  adduct 
homopolymerization  via  propenyl  groups,  and  dehydration/ 
etherification  reaction)  leading  to  a  much  more  complex  curing 
pathway.  '1’  ’"  The  framework  for  the  curing  process  involving 
all  five  reactions  is  currently  being  formulated  using  atomistic 
molecular  dynamics  simulations  and  would  be  reported  in 
future  studies.  Nevertheless,  based  on  our  previous  studies  on 
interface  thermal  conductance  across  cured  epoxy  matrix  and 
carbon  nanotubes,  we  foresee  that  the  curing  process  is 
expected  to  enhance  the  thermal  conductance  up  to  20%. “4 
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4.  CONCLUSIONS 

In  a  single  comprehensive  investigation,  this  study  explores  a 
series  of  carbon-fiber  bulk  and  interface  parameters  in  order  to 
gauge  their  importance  toward  engineering  thermal  transport 
across  carbon-fiber/matrix  interfaces.  Based  on  atomistic 
molecular  dynamics  simulations,  it  is  suggested  that  the  long- 
range  graphitic  nature  of  carbon-fibers  near  the  interface,  a 
higher  number  density  of  surface  atoms,  a  larger  degree  of 
functionalization,  molecular  surface  roughness,  and  incorpo¬ 
ration  of  high-conductivity  fillers  (such  as  CNTs)  can 
significantly  augment  the  interface  thermal  conductance  across 
the  interface.  On  the  contrary,  defects  within  the  carbon  fiber, 
low  number  density  of  interacting  atoms,  and  less  effective  area 
of  interaction  due  to  voids  (not  studied  here)  can  lead  to  a 
noticeable  decrease  in  thermal  conductance.  In  summary,  the 
presented  results  provide  key  insights  into  realistic  bounds 
(degree  of  variation)  to  the  interface  thermal  conductance 
values  at  fiber/matrix  interfaces  as  a  function  of  many  different 
surface-carbon  morphologies.  In  addition,  we  believe  that  this 
parametric  study  incorporating  various  key  parameters  has 
broader  ramifications  on  modulating  thermal  conductance,  not 
only  across  carbon-fiber/matrix  interfaces  but  also  in  different 
types  of  thermal  interface  materials  as  well  as  different 
nanoscopic  devices  where  thermal  dissipation  plays  a  crucial 
role  in  device  performance. 
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SI:  Generation  of  carbon  surface  models  with  different  graphitic  orientations 


Figure  SI:  Carbon  surface  models  with  (a)  0°  graphitic  orientation;  (b)  90°  graphitic  orientation;  and  (c)  33° 
graphitic  orientation;  (d)  close-up  view  highlighting  periodic  nature  of  the  titled  graphitic  layer  as  wrapped 
red  lines;  (e  and  f)  mathematical  constraints  applicable  for  generating  periodic  structures  with  non-0°  and 
non-90°  graphitic  orientations;  and  (g)  final  graphitic  tilt  angles  generated  from  the  mathematical  equations. 

To  investigate  the  effect  of  graphitic  orientation  on  the  interface  thermal  conductance,  several  models 
were  created  having  different  graphitic  orientations  with  respect  to  the  interface  plane  (horizontal  plane  in 
Figure  SI).  It  is  relatively  easier  to  create  parallel  and  perpendicular  orientations  of  the  graphitic  periodic 
crystal,  as  shown  in  Figures  Sla  and  Sib  using  graphitic  unit  cell.  However,  creating  a  ‘periodic  struc¬ 
ture’  with  any  other  graphitic  orientation  (as  shown  in  Figure  Sic)  is  not  trivial  and  certain  constraints 
need  to  be  satisfied  towards  their  generation.  These  constraints  are  shown  in  Figures  Sld-Slf  and  are 
discussed  below. 

Graphite  is  periodic  along  zigzag  and  armchair  directions  with  periodic  lengths  of  2.46  and  4.26  A,  re¬ 
spectively.  Towards  the  generation  of  graphitic  orientations  (as  shown  in  Figure  Sic),  we  chose  to  have 
armchair  direction  along  the  Y-axis  (going  into  the  plane)  while  the  zigzag  direction  makes  certain  angle 
with  the  X-axis  (horizontal  axis  in  the  plane).  So,  Figure  Sla  has  the  orientation  angle  of  0°  while  Figure 
Sib  has  the  orientation  angle  of  90°  with  zigzag  direction  being  horizontal  and  vertical  (in  plane  of  pa¬ 
per),  respectively.  We  call  them  C()  and  C90  systems  respectively. 

First,  along  y-diicction,  12  periodic  cells  (along  armchair)  were  used,  rendering  the  y-dimension  of  ~  5.1 
nm.  In  addition,  x-dimension  length  was  set  to  be  5  nm  approximately. 
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One  constraint  of  creating  non-0°  and  non-90°  graphitic  orientations  is  that  since  individual  graphene  lay¬ 
ers  are  periodic,  one  graphene  layer  appears  as  several  graphene  layers  (separated  by  a  certain  distance) 
when  wrapped  across  the  periodic  box.  Such  wrapping  behavior  is  shown  in  Figure  Sid  (red  diagonal 
lines  which  are  indicative  of  a  single  graphene  layer  wrapped  around  the  periodic  boundary).  Another 
constraint  that  one  might  have  is  the  AB  stacking  geometry  of  graphene  layers  within  graphite.  When  we 
take  both  constraints  into  account,  it  can  be  realized  that  eventual  graphitic  orientation  angles  will  be  set 
by  certain  trigonometric  equations  and  it  is  not  possible  to  have  any  angle  if  we  constrain  X-dimension  to 
be  ~5  nm.  These  trigonometric  equations  are  shown  in  Figure  Sle. 

Looking  at  the  triangle  in  Figure  Sle  and  setting  its  horizontal  length  to  be  5  nm,  angle  0  will  be  decided 
by  how  many  different  graphitic  layers  we  have  in  between  before  the  same  layer  repeats  itself  (i.e.,  lay¬ 
ers  between  red  lines).  W,  width  perpendicular  to  the  graphene  plane  will  be  given  by  0.68/7  for  AB  type 
of  stacking,  where  n  is  number  of  AB  repeat  units  perpendicular  to  the  graphitic  plane.  Using  these  two 
constraints,  we  can  find  out  angle  0. 

However,  we  have  to  fulfill  one  more  constraint.  The  length  L  which  is  determined  by  5/cos(0)  cannot  be 
any  length.  This  length  has  to  be  a  multiple  of  periodicity  (2.46  A)  along  zigzag  direction  (shown  as  /.). 
To  have  all  constraints  valid,  we  wrote  a  script  to  come  up  with  different  graphitic  orientations.  These 
orientations  would  have  well-defined  angles  based  on  the  value  of  n  and  horizontal  distance  (5  nm  in  our 
case)  which  we  put  into  the  script.  These  orientation  angles  are  shown  in  the  table  form  in  Figure  Slg. 

Using  these  constraints,  7  additional  graphitic  models  (in  addition  to  0  and  90)  were  created  for  n  ranging 
from  1  to  7.  In  all  cases  towards  studying  the  effect  of  graphitic  orientation,  the  thickness  of  the  carbon 
surface  model  was  fixed  to  be  10  nm.  Systems  were  equilibrated  as  discussed  in  the  main  text. 


S-3 

36 

Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 


S2:  Generation  of  interface  models  with  different  functionalization 


Figure  S2:  (a)  BMI  monomer:  4,4’-bismaleimidodiphenyI  methane  (BMPM);  (b)  BMI  monomer:  0,0’  diallyl 
bisphenol  A  (DABPA);  (c)  combined  starting  system  of  C90  carbon  surface  model  and  equilibrated  BMI  mon¬ 
omeric  matrix;  (d)  close-up  view  of  starting  confirmation  of  2  nm  region  close  to  the  interface;  (e)  close-up 
view  of  final  confirmation  of  2  nm  region  after  dragging  simulations;  (f)  top  view  of  the  interface  ‘O’  atoms; 
(g)  perspective  view  of  interface  ‘O’  atoms. 

Figure  S2  shows  the  first  phase  of  functionalization  model  that  we  used  for  functionalizing  the  surface  of 
Cg0  graphitic  model.  Figures  S2a  and  S2b  show  the  monomeric  units  for  the  BMI  matrix  that  we  have 
used  for  modeling  the  matrix.  First,  the  equilibrated  matrix  is  combined  with  the  C90  graphitic  model  as 
shown  in  Figure  S2c  using  an  in-house  script  in  LAMMPS  format.  Moreover,  ‘O’  atoms  from  hydroxyl 
groups  of  monomer  DAPBA  (red  atoms  in  Figure  S2b)  in  lower  2  nm  region  are  selected  and  dragged 
towards  the  interface  for  subsequent  bonding.  These  selected  atoms  are  shown  as  red  spheres  in  Figure 
S2c-S2g.  The  close-up  view  of  a  2  nm  region  in  the  initial  state  is  shown  in  Figure  S2d.  A  MD  simula¬ 
tion  is  subsequently  run  where  a  force  of  1.0  Kcal/mole -Angstrom  (LAMMPS  units)  is  applied  on  each 
atom  to  drag  them  near  the  interface  for  200  ps  with  a  time-step  of  0.5  fs.  Figure  S2e  shows  the  final 
structure  of  the  2  nm  region  after  the  drag  simulations.  Figures  S2f  and  S2g,  respectively  show  the  top 
and  perspective  view  of  the  ‘O’  atoms  which  are  very  close  to  the  surface. 

Once  the  ‘O’  atoms  are  near  the  surface,  the  nearest  carbon  atoms  (from  C90  model)  for  each  oxygen  atom 
is  identified  for  potential  reaction  (up  to  13  atoms).  These  atoms  are  identified  in  a  sequence  for  reaction 
such  that  functionalizing  atoms  are  as  far  away  as  possible.  This  helps  in  maximizing  the  distribution  of 
the  functionalization  groups.  Towards  this,  an  in-house  script  was  written  to  identify  the  sequence  of  at¬ 
oms  for  reactions  which  maximizes  inter-particle  distance  for  cases  when  we  have  only  few  functional- 
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ized  atoms.  This  sequence  is  shown  in  Figure  S3  in  part  for  1,  4,  7,  10  and  13  functionalized  carbon  at¬ 
oms  leading  to  functionalization  from  0.05  nm 2  to  0.55  nm  2.  Please  note  that  the  surface  is  in  fact  period¬ 
ic  in  lateral  X-Y  dimensions. 


Carbon  edge  atoms  identified  for  reaction  with  increasing  functionalization 


Figure  S3:  Schematic  of  5  images  showing  the  progression  of  selected  carbon  atoms  for  functionalization  as 
we  increase  functionalization  surface  density,  calculated  using  a  script  which  maximizes  distance  between  the 
bonded  carbon  atoms. 


Figure  S4:  (Right)  Initial  structure  of  the  interface  before  bonding;  (Left)  Final  structure  of  the  interface  af¬ 
ter  bonding  for  13  atom  functionalization  (functionalization  surface  density  of  0.55  nm'2). 

Once  the  atoms  are  identified,  we  create  a  bond  between  O  (red)  and  C  (blue)  atoms  and  update  the  topol¬ 
ogy  information  of  the  LAMMPS  datafile  (new  bonds,  angles,  dihedrals,  impropers)  and  minimize  the 
newly  bonded  structure.  The  initial  and  final  geometry  is  shown  in  Figure  S4. 

Subsequent  equilibration  simulations  are  run  as  discussed  in  the  main  text,  similar  to  other  cases. 
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S3:  Generation  of  interface  models  with  different  surface  roughness 


Sinusoidal  Roughness  Parameters: 
Amplitude  (A),  Wavelength  (W) 


(c) 


Figure  S5:  (a)  Starting  C90  carbon  surface  model;  (b)  Schematic  of  how  the  region  is  selected  to  delete  atoms 
for  creating  roughness  of  amplitude  A  and  wavelength  W;  (c)  Schematic  representation  of  created  rough  C90 
carbon  surface  models  in  side  and  perspective  views  (top  and  bottom  images,  respectively).  Bottom  images 
also  include  periodic  images  of  the  created  rough  surfaces  (3X3  periodic  boxes  along  X  and  Y  direction)  in 
initial  configuration 


Figure  S5  shows  the  steps  for  creating  roughness  in  the  modeled  carbon  surface  models.  For  these  set  of 
simulations,  we  used  a  C90  model  as  shown  in  Figure  S5a.  The  reasoning  for  using  C90  is  discussed  later. 
In  order  to  create  roughness  features  of  amplitude  A  and  wavelength  W  (as  shown  in  Figure  S5b),  we 
create  a  fictitious  plane  (parallel  to  X-Y  plane)  A  A  (amplitude)  underneath  the  surface  as  shown  by  the 
horizontal  double  sided  arrow.  Then,  we  create  10000  (100  bins  X  100  bins)  rectangular  slabs  of  different 
heights  (based  on  2D  sinusoidal  curve)  using  the  following  script  in  LAMMPS  (using  lammps  region 
module)  as  shown  in  partially  shaded  orange  region  in  Figure  S5b. 


in . cut 


variable 

variable 

variable 

variable 

variable 


zmax  equal  bound ( all , zmax) 
center  equal  ${ zmax} -0 . 5*$ {amp} 
lowx  equal  xlo 
xd  equal  lx 

binlength  equal  $ { xd} / ( $ {bin } *  1 . 0 ) 


label 

variable 

variable 


loop 

a  loop  ${bin} 

tempx  equal  ${ lowx }+((${ a } *1 . 0-0 . 5 )*${ binlength} ) 
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variable 

variable 

variable 

variable 

print 

if  "${zmax}  > 

"region 

"delete^atoms 

"region 

variable 

variable 

variable 

variable 

variable 

next 

j  ump 


xl  equal  $ { lowx } + ( ( $ {a } *  1 . 0-1 . 0) *$ {binlength } ) 
x2  equal  ${ lowx }+((${ a }* 1 . 0) *$ {binlength } ) 
xx  equal  2* 4*atan ( 1 . 0 ) *$ { tempx } / $ { xd} 
zl  equal  ${ center } +0 . 5*$ {amp } *sin ($ {period} *${ xx} ) 

"$a  ${zl}  $ { zmax } " 

${zl}"  then  & 

temp  block  ${xl}  ${x2}  EDGE  EDGE  ${zl}  ${zmax}  units  box"  & 

region  temp"  & 

temp  delete" 

tempx  delete 

xl  delete 

x2  delete 

xx  delete 

zl  delete 

a 

in. cut  loop 


33°  GRAPHITIC  ORIENTATION  MODEL 


Along  Y  direction  Along  X  direction 


Perspective  View  Top  View 


0°  GRAPHITIC  ORIENTATION  MODEL 


Figure  S6:  Right:  C33  carbon  surface  model  showing  roughness  features  as  viewed  along  x-  and  y-direction. 
The  bottom  images  show  the  perspective  and  top  view  of  the  rough  C33  interface;  Left:  C0  carbon  surface 
model  as  viewed  along  j -direction  (into  the  plane). 

Here,  the  top  z-coordinate  (zmax)  of  the  individual  slab  is  constant  while  the  bottom  z-coordinate  (zl)  var¬ 
ies  based  on  the  equations  above.  Moreover,  slab  bins  along  x-  and  y-directions  are  determined  by  overall 
box  length  along  the  respective  directions,  divided  by  100.  Once  we  have  all  regions  defined,  we  delete 
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the  atoms  in  these  regions,  resulting  in  sinoduidal  roughness  features  on  the  surface  as  shown  in  Figure 
S5c  for  different  amplitude  and  wavelengths.  Thereafter,  we  further  delete  a  few  more  atoms  which  have 
only  1  atom  connected  resulting  from  the  previously  deleted  atoms.  After  that,  we  update  the  topology 
information  (removing  corresponding  bonds,  angles,  dihedrals  and  improper  angles)  to  generate  a  final 
LAMMPS  data  file  for  rough  C90  carbon  surface  models,  followed  by  their  subsequent  equilibration  as 
discussed  in  the  main  text. 

Figure  S6  depicts  the  reason  for  using  only  C90  model  for  roughness  simulations.  The  figure  shows  2  car¬ 
bon  surface  models  (C33  and  Co)  with  created  roughness  features.  As  can  be  seen  from  the  figure,  the  gen¬ 
erated  roughness  features  also  lead  to  generation  of  small  un-connected  fragments  (top  region).  It  is  ex¬ 
pected  that  these  un-connected  fragments  would  diffuse  into  the  matrix  during  the  combined  system’s 
equilibration,  thus  defying  the  sole  purpose  of  studying  the  effect  of  surface  roughness  on  interface  ther¬ 
mal  transport.  Hence,  for  investigating  rough  systems,  other  orientations  except  C90  were  not  considered. 
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S4:  Generation  of  interface  models  with  CNTs  incorporated  at  the  interface 


Figure  S7:  (a)  Schematic  of  (10,10)  SWCNT  fuzed  into  C90  graphitic  surface  carbon  model;  (b)  close-up  view 
of  the  interface;  (c)  perspective  view  of  the  close-up  interface;  (d)  perspective  view  of  the  capped  (10,10) 
SWCNT;  (e)  close-up  view  of  the  SWCNT  cap  in  equilibrium  state;  (f)  schematic  of  BMI  matrix,  showing  a 
cylindrical  hole  in  the  BMI  matrix  for  SWCNT  inclusion;  (g)  starting  geometry  of  the  combined  matrix  and 
fuzed-SWCNT  carbon  surface  model;  and  (h)  snapshot  of  the  combined  system  after  100  ps  NPT  simulation 
showcasing  that  most  of  the  empty  space  is  disappeared. 

In  order  to  create  a  fused  junction  of  carbon  nanotube  (CNT)  with  graphitic  carbon  surface  model,  the  C,m 
surface  model  was  chosen.  This  model  was  primarily  chosen  for  the  reason  that  CNTs  could  be  fused  per¬ 
pendicular  to  the  interface  with  the  most  connectivity  (most  CNT  atoms  connected  to  the  C90  surface  at¬ 
oms).  One  such  overall  junction  is  shown  in  Figure  S7a  (junctions  with  different  CNT  lengths  were  cre¬ 
ated).  The  close-up  version  of  the  interface  is  shown  in  Figure  S7b.  This  sub-figure  also  depicts  why 
(10,10)  CNT  was  chosen  to  fuse  with  the  C90  graphitic  model.  As  seen  from  the  figure,  the  diameter  of 
(10,10)  is  13.6  A,  which  is  very  similar  to  the  distance  between  4  graphitic  layers.  This  matching  resulted 
in  12  of  the  20  atoms  getting  connected  to  4  graphitic  layers.  Also,  the  symmetry  of  the  nanotube  (along 
C2V  axis)  as  well  as  its  arm-chair  nature  resulted  in  perfect  connectivity  at  both  ends  of  the  CNT.  This 
smooth  transition  of  connectivity  is  shown  in  Figure  S7c,  where  6  benzene  rings  (3  on  each  side)  are 
smoothly  transitioned  from  the  graphite  layers  to  CNT.  Once  the  bonds  were  created,  the  new  topology 
information  (new  bond  angles,  new  dihedrals,  new  improper  angles)  was  incorporated  in  the  updated 
LAMMPS  data  file  using  an  in-house  script. 

During  the  preparation  it  also  occurred  to  us  that  because  of  the  large  diameter  of  fused  CNTs  (>  lnm), 
BMI  matrix  molecules  might  get  trapped  inside  the  CNTs  because  of  its  open  nature  on  one  end  (end 
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away  from  the  graphitic  interface).  In  order  to  circumvent  the  problem,  (10,10)  capped  SWCNTs  (on  one 
end)  were  created  and  fused  at  the  open  end  using  the  methodology  discussed  above.  Figure  S7d  shows 
the  graphitic  carbon  surface  model  with  a  fused-capped-SWCNT  while  Figure  S7e  shows  the  close-up  of 
the  equilibrated  nanotube  cap. 

In  order  to  combine  the  CNT  incorporated  carbon  surface  C90  model  into  the  BMI  matrix,  a  20  A  diameter 
hole  was  created  within  the  BMI  matrix  along  the  z-direction  as  shown  in  Figure  S7f.  This  was  done  us¬ 
ing  “fix  indent”  module  of  LAMMPS  MD  software.  The  hole  was  radially  grown  from  a  radius  of  0  to  10 
A  over  100  ps  MD  run  to  avoid  sudden  changes  in  the  geometry  of  the  BMI  matrix  molecules.  Once  the 
hole  was  created,  both  systems  were  combined  in  a  single  LAMMPS  data  file  using  an  in-house  script  as 
depicted  in  Figure  S7g  and  were  equilibrated  for  several  nanoseconds  (Figure  S7h). 

Using  this  methodology,  4  different  systems  with  CNT  length  ranging  from  3  nm  to  9  nm  were  created  to 
understand  the  effect  of  CNT  length  (and  thus  effective  surface  area)  on  interfacial  thermal  transport 
across  BMI-matrix  carbon-fiber  interface. 
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S5:  Estimation  of  temperature  gradient  across  the  interface  using  steady-state  tem¬ 
perature  profile. 
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Figure  S8:  (Main)  Steady-state  temperature  profile  of  a  representative  case  showcasing  the  temperature  drop 
across  the  interface  of  carbon-fiber  and  matrix;  (Inset)  Steady-state  temperature  profile  of  the  carbon-fiber. 

Figure  S8  shows  the  steady-state  temperature  profile  as  well  as  the  temperature  drop  across  the  interface. 
The  inset  also  shows  that  the  monotonic  temperature  drop  (amid  fluctuations)  across  the  carbon  fiber  is 
very  small  since  the  thermal  conductivity  of  the  graphitic  carbon  fiber  is  significantly  higher  than  that  of 
matrix.  In  order  to  calculate  the  temperature  drop,  the  data  from  both  side  of  the  interface  was  fitted  to  a 
linear  fit  as  shown  in  the  figure.  Then,  based  on  the  fitting,  AT  at  the  interface  was  calculated  as  depicted 
in  the  figure  for  thermal  conductance  estimation. 
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S6:  Estimation  of  heat  flux  using  LAMMPS  output  data 


Figure  S9:  Plot  of  cumulative  energy  added  (blue)  and  taken  out  (green)  from  hot  and  cold  thermostats,  re¬ 
spectively,  as  a  function  of  simulation  time  for  a  representative  case.  The  slopes  are  linear  fit  to  data. 

Figure  S9  shows  the  LAMMPS  output  cumulative  energy  as  a  function  of  simulation  time.  The  heat  flux 
is  calculated  from  the  slope  by  fitting  the  data  to  a  linear  fit,  followed  by  its  division  by  the  perpendicular 
surface  area  (along  X-Y  plane),  thus  giving  the  units  of  kcal/mol/A2/ps,  which  was  subsequently  convert¬ 
ed  to  W/m2. 
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S7:  Effect  of  using  various  interfaces  for  interface  thermal  conductance  calculations 

for  rough  interfaces 


Figure  S10:  (Top)  Figure  showing  the  overlap  region  of  the  rough  interface  in  terms  of  its  temperature  calcu¬ 
lations.  (Bottom)  Plots  of  interface  thermal  conductance  as  a  function  of  RMS  roughness  for  different  period¬ 
ic  roughness  values,  calculated  at  directed  interfaces. 

Figure  S10  clearly  depicts  the  ambiguity  of  calculating  ‘interface’  thermal  conductance  across  rough  in¬ 
terfaces.  As  seen  from  the  figure,  the  calculations  provide  very  different  quantitative  values  as  well  as 
qualitative  trends  depending  upon  the  chosen  interface.  The  reason  being  that  although,  the  heat  flux  and 
cross-sectional  area  are  constant,  AT  is  noticeably  different  at  selected  interfaces.  In  such  cases,  it  be¬ 
comes  ambiguous  as  to  which  interface  should  be  chosen.  To  avoid  that  ambiguity,  we  calculated  ‘inter¬ 
phase’  thermal  conductance  between  points  A  and  B  as  discussed  in  the  main  text. 
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S8:  Visual  differences  between  C0  and  C$-C90  graphitic  orientation  cases  with  respect 

to  applied  thermostats 


c 

o 


Figure  SI  1 :  (Right)  C0  graphitic  orientation  model;  (Left)  C4 1  graphitic  orientation  model.  Here,  3  periodic 
images  (with  periodic  boxes)  of  the  graphitic  layers  are  also  shown  to  further  clarify  the  periodically  wrapped 
graphene  connectivity. 

Figure  Sll  visually  shows  the  difference  of  how  heat  is  transported  from  hot  thermostats  to  graphitic  in¬ 
terface  region  in  C0  vs  C4[  graphitic  orientation  cases.  For  the  C0  model,  heat  propagates  along  the  out-of¬ 
plane  direction  of  graphitic  layers  before  reaching  the  interface.  On  the  other  hand,  for  C«  model,  heat 
propagates  parallel  to  the  graphitic  plane  and  reaches  the  interface  due  to  the  huge  thermal  anisotropy  of 
graphite  (thermal  conductivity  along  graphite  plane  is  >100  times  higher  than  out-of-plane  direction).  We 
should  also  point  out  that  the  effective  length  of  heat  propagation  here  is  significantly  longer  than  1 0  nm 
and  is  determined  by  the  orientation  angle. 
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ABSTRACT 

Very  rapid  heating  of  composite  material  leads  to  matrix  thennal  decomposition, 
amorphous  char  formation,  ablation,  etc.  and  results  in  the  variability  of  local  thermo¬ 
physical  properties.  It  is  desired  to  delay  the  onset  of  matrix  degradation  by 
incorporating/tailoring  the  fiber-matrix  interface  for  efficient  thennal  energy  transfer 
and  helping  in  faster  dispersion/distribution  of  supplied  thennal  energy  and  thus, 
necessitates  the  understanding  of  the  interface  as  well  as  interfacial  thennal  transfer 
across  the  fiber/matrix  interface.  Within  the  context  of  modeling  laser  inadiation 
based  rapid  heating  of  nanocomposites,  i.e.,  modeling  of  thennal  energy  transfer 
across  the  composite’s  interfaces  and  subsequent  polymeric  ablation,  this  study 
focuses  on  building  a  model  and  investigating  physical  and  thennal  properties  of  high- 
temperature  imide-based  resin  (BMI-5292)  matrix  using  molecular  dynamics 
simulations.  A  molecular  model  of  BMI-5292  based  on  AMBER  force  field  and 
GAUSSIAN-derived  RESP  charges  is  developed  and  validated  against  several 
experimentally  measured  thenno-physical  properties  such  as  density,  glass  transition, 
thennal  conductivity,  etc.  In  parallel,  to  mimic  the  carbon  fiber  surface,  a  model  of 
amorphous  carbon  is  developed  using  the  Tersoff  (bond-order)  potential. 
Subsequently,  using  non-equilibrium  molecular  dynamics  simulations,  interfacial 
thennal  energy  exchange  is  investigated  in  terms  of  thennal  conductance  at  carbon 
fiber-matrix  interfaces  for  different  degrees  of  crystallinity  within  the  carbonaceous 
region.  The  obtained  thenno-physical  results  are  found  to  be  in  good  agreement  with 
available  literature,  thus  providing  an  accurate  and  reliable  foundation  of  the 
molecular  force-field  to  tackle  more  complex  phenomenon  (such  as  polymeric 
ablation)  from  molecular  dynamics  simulations  with  respect  to  high-temperature 
resins. 
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INTRODUCTION 


Today,  laser  technology  is  being  used  in  a  wide  variety  of  fields  including 
communications,  industrial  and  enviromnental  applications,  medicine,  and  materials’ 
research  and  development  (R&D).  Some  of  the  key  areas  in  materials’  R&D  include 
mass  spectral  analysis  using  matrix-assisted  laser  desorption/ionization  (MALDI)  and 
fabrication  of  micro-/nano-pattemed  surfaces  and  thin  films  using  pulsed  or 
continuous  lasers  (such  as  micro-fabrication  of  electronic  devices  and  micro- 
patterning  of  polymeric  composites).  Our  specific  interest  lies  in  how  the  different 
lasers  interact  with  polymeric  composite  materials;  what  parameters  are  of  important 
significance  which  governs  their  mutual  interaction;  and  how  such  interaction  leads  to 
matrix  thennal  decomposition,  amorphous  char  fonnation,  and  polymeric  ablation. 

The  basic  mechanism  of  laser  interaction  with  materials  incurs  many  non¬ 
equilibrium  photo-processes  (thermal,  mechanical,  chemical,  physical)  caused  by  fast 
deposition  of  laser  energy  and  its  dissipation  in  various  forms  [1,  2].  Because  of  the 
rapid  heating  rate,  these  non-equilibrium  processes  cause  spatial  variability  in  thermo- 
physical  properties  of  the  matrix  of  interest  such  as  local  heating,  local  matrix 
degradation  due  to  sudden  local  increase  in  temperature,  pressure  wave  fonnation,  etc. 
[1,  2].  It  is  often  desired  to  delay  the  onset  of  ablation  and  keep  the  composite  intact  by 
tailoring  the  interfaces  of  matrix  and  fillers  (such  as  carbon  fibers,  carbon  nanotubes, 
etc.)  by  modulating  interactions  between  the  two  constituents.  To  better  understand 
how  lasers  interact  with  composites,  one  needs  to  investigate  the  correlation  between 
experimental  laser  parameters,  i.e.,  type  of  laser  (pulsed  or  continuous),  its  fluence  and 
incident  angle,  penetration  depth,  and  composite  materials  parameters,  i.e.,  molecular 
chemistry,  cohesive  energy,  thenno-physical  properties,  mechanical  properties, 
matrix-filler  interfacial  interactions,  etc. 

Figure  1  shows  the  occurrence  of  different  physical  and  chemical  events  after  laser 
absorption  along  with  their  respective  timeframe  [1].  As  these  events  span  over  several 
orders  of  magnitude  in  time,  it  is  impossible  to  provide  a  consistent  analytical 
description  of  all  processes  and  address  all  relevant  physics  and  chemistry  within  a 
single  computational  framework.  It  can  be  seen  from  the  figure  that  a  multitude  of 
interesting  molecular  events  occur  at  the  picosecond  (ps)  to  nanosecond  (ns)  time 
scales,  investigations  and  understanding  of  which  is  best  suited  using  molecular 
dynamics  (MD)  simulations  [3]  as  it  operates  at  these  timescales.  In  these  simulations, 
only  the  details  of  microscopic  interactions  need  to  be  specified  (often  extracted  using 
quantum  chemical  calculations)  and  no  further  assumptions  need  to  be  made  in  order 
to  investigate  the  equilibrium  and  non-equilibrium  phenomenon  happening  at  these 
timescales. 

With  our  specific  interest  in  Matrimid™  BMI-5292  matrix  based  carbon-fiber 
composites  and  to  understand  how  laser  absorption  leads  to  the  heat  dissipation, 
matrix  degradation  and  the  eventual  ablation  within  our  system  of  interest,  this  study 
provides  an  initial  framework  along  the  specified  questions.  Specifically,  in  this  study, 
we  built  a  full-atomistic  model  of  BMI-5292  resin  monomer  mixtures  and  calculate 
their  thenno-physical  properties  (density,  glass  transition  temperature,  thennal 
conductivity,  coefficient  of  linear  thennal  expansion)  using  MD  simulations  and 
compare  the  results  with  reported  experimental  values.  Furthennore,  we  investigate 
the  thennal  energy  exchange  across  the  carbon  fiber  (amorphous) — BMI-matrix 
interface.  We  believe  that  this  initial  investigation  provides  a  foundation  for  the 
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validity  and  accuracy  of  the  employed  MD  force-field  and  can  be  built  upon  further  to 
investigate  more  complex  issues  of  interests  such  as  polymeric  ablation  with 
molecular  precision. 
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Figure  l.A  schematic  flow  of  events  that  occur  after  absorption  of  laser  energy  along  with  their 
respective  timeframes.  The  partially  transparent  shaded  region  (~1  ps  -10s  of  ns)  is  the  timeframe  that  is 

best  suited  for  MD  simulations. 


SIMULATION  METHODOLOGY 
Simulation  Setup  and  Methodology:  BMI  Matrix 

The  BMI  resin  (Matrimid™  BMI-5292)  is  a  copolymer  composed  of  4,4’- 
bismaleimidodiphenyl  methane  (BMPM)  and  0,0’  diallyl  bisphenol  A  (DABPA) 
monomers  in  a  1:1  molar  ratio  [4],  The  schematic  representations  of  these  monomers 
are  shown  in  Figure  2.  As  a  part  of  simulation  set-up,  these  monomers  were 
constructed  in  Materials  Studio®  Visualizer  [5].  Several  force-fields  (CVFF  [6],  PCFF 
[7],  AMBER  [8,  9])  were  tested  for  successful  generation  of  bonded  and  non-bonded 
parameters.  However,  only  the  general  amber  force  field  (GAFF)  was  able  to  generate 
all  required  parameter  sets  among  the  tested  force  fields  [8],  Hence,  the  GAFF  force- 
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field  was  employed  for  testing  its  accuracy  and  to  model  BMI  matrix  monomers  for 
this  study. 

It  is  often  reported  that  when  using  the  GAFF  force-field,  restrained  electrostatic 
potential  (RESP)  charges  are  best  suited  to  model  organic  and  bio-molecules  [8,  9]. 
For  our  simulations,  RESP  charge  estimations  and  fittings  were  perfonned  using 
GAUSSIAN09  software  [10]  and  RESP  code  within  the  framework  of  the  R.E.D. 
program  [11].  Further  details  on  the  charge  estimation  will  be  provided  later  in  the 
results  and  discussion  section.  Once  all  the  required  parameterization  was  done,  the 
GAFF  force-field  and  the  coordinates  of  the  BMI  matrix  components  were  imported 
into  LAMMPS  molecular  dynamics  package  to  carry  out  molecular  dynamics  (MD) 
simulations  [12],  A  cutoff  of  12  A  was  used  for  van-der  Waals  interactions  and  short 
range  electrostatics,  while  long  range  electrostatics  was  treated  using  the  PPPM 
methodology  [13]  as  implemented  in  LAMMPS.  Visual  molecular  dynamics  (VMD) 
software  was  employed  to  visualize  the  trajectories  and  create  representative 
schematic  visualizations  that  are  reported  in  this  study  [14]. 


o  o 

n  ii 


4,4-bismaleimidodiphenyl  methane 


0,0-diallyl  bisphenol  A 


Figure  2.  Molecular  representation  of  the  two  monomer  components  of  the  Matrimid™  BMI-5292  resin. 
Color  Scheme:  cyan-carbon;  blue-nitrogen;  red-oxygen;  silver-hydrogen. 


First,  a  system  of  -9500  atoms  was  created  by  replicating  both  BMI  monomers 
108  times  (total)  along  x-,  y-,  and  z-directions  as  shown  in  Figure  3(a).  Then,  MD 
simulations  were  carried  out  in  NVT  (canonical)  ensemble  to  equilibrate  the 
temperature,  followed  by  NPT  (isothermal-isobaric)  ensemble  to  equilibrate  the 
density  (and  pressure).  Periodic  boundary  conditions  were  employed  in  all  3 
orthogonal  directions  to  mimic  bulk-like  behavior.  While  the  successful  temperature 
equilibration  (NVT  simulations)  was  achieved  within  200  ps,  the  equilibration  of 
density  (NPT  simulations)  required  significantly  longer  simulations.  In  total,  -10  ns  of 
NPT  simulations  with  independent  barostats  along  the  3  directions  were  performed  for 
successful  equilibration  of  density.  Afterwards,  a  relatively  shorter  run  of  200  ps  was 
perfonned  in  NVE  (micro-canonical)  ensemble  to  confirm  the  energy  conservation  of 
the  system.  A  timestep  of  1.0  femtosecond  (fs)  was  used  for  these  set  of  simulations 
during  the  equilibrium  phase.  The  representative  equilibrated  snapshot  is  shown  in 
Figure  3(b).  After  successful  equilibration,  the  system  was  subjected  to  a  series  of 
additional  simulations  for  estimation  of  glass  transition  temperature  (Tg),  density 
variations  with  temperature,  coefficient  of  linear  thermal  expansion  (CLTE)  and 
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thermal  conductivity  (A).  For  the  sake  of  continuity  of  discussion,  the  methodology  for 
the  specific  simulations  will  be  summarized  in  appropriate  sections  later  on. 


(a)  (b) 

Figure  3.  Schematic  representation  of  the  (a)  initial  and  (b)  equilibrated  geometry  of  the  BMI-5292 
component  mixture  (-9500  atoms)  as  discussed  in  the  text.  Visually,  the  red-bluish  and  green-blackish 
molecules  represent  BMPM  and  DABPA,  respectively. 


Simulation  Setup  and  Methodology:  Carbon  Fiber  Surface 

Initially,  a  diamond  lattice  of  dimensions  -5x5x5  mn  was  generated  using 
Materials  Studio®  Visualizer.  The  structure  was  then  rapidly  heated  to  12000  K  in 
NVT  ensemble  using  the  Tersoff  potential  [15]  as  implemented  in  LAMMPS  MD 
package.  Periodic  boundary  conditions  were  used  in  all  3  directions  at  this  stage.  Due 
to  the  short-range  nature  of  the  Tersoff  potential,  the  large  temperature  fluctuations  at 
elevated  temperatures  led  to  randomization  of  the  lattice,  and  thus,  the  generation  of 
the  amorphous  structure.  This  amorphous  structure  was  then  allowed  to  cool  down  to 
room  temperature  at  three  different  cooling  rates  (lK/ps,  1  OK/ps  and  lOOK/ps)  in  NPT 
ensemble  with  independent  barostats  in  all  3  periodic  directions.  During  cooling, 
nano-crystalline  diamond  domains  were  fonned  for  the  low  cooling  rate  (lK/ps), 
while  a  fully  amorphous  structure  evolved  for  higher  cooling  rates.  Afterwards,  the 
cooled  structures  were  further  relaxed  under  NPT  and  NVE  ensembles  for  -500  ps  at 
room  temperature  for  density  and  energy  conservation,  respectively.  The  flow  of  the 
carbonaceous  system  simulations  is  better  represented  in  Figure  4,  along  with 
representative  snapshots  of  annealed  and  equilibrated  bulk  structures.  Here,  we  should 
point  out  that  since  the  Tersoff  potential  does  not  have  a  dispersion  interaction 
component  to  it,  graphitic  order  (enhanced  sp2  order  and  graphene  layers)  was  not 
observed  in  nano-crystalline  domains.  Once  the  fully  relaxed  amorphous  bulk-like 
carbon  structures  were  realized,  two  surfaces  (for  each  case  of  annealed  carbon 
structure)  were  generated  by  increasing  the  z-dimension  of  the  box.  The  structure  was 
further  relaxed  in  NPT  ensemble  at  room  temperature  to  relax  the  stresses  of  the 
newly-generated  surfaces.  A  time  step  of  0.5  fs  was  employed  throughout  the 
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equilibrium  stage.  We  should  point  out  that  another  bond-breaking  potential,  namely 
AIREBO  [16],  was  also  tried  to  generate  equilibrated  structures  of  the  matrix. 
However,  none  the  simulations  were  successfully  finished  and  were  crashed  during  the 
initial  heating  stage  of  the  process. 


Diamond  Lattice 


Heated  amorphous  Structure 


Increasing  annealing  rate 


r 

Final  annealed  structures,  annealed  at  different  rates 


Figure  4.  (Top  Left)  Initial  diamond  lattice;  (Top  Right)  Heated  amorphous  structure  @6000  K; 
(Bottom)  Equilibrated  bulk  structures,  obtained  by  cooling  down  the  top-right  structure  with  different 
rates.  It  is  clear  the  fast  cooling  leads  to  amorphous  structure  while  slower  cooling  lead  to  the  formation 

of  nano-crystalline  diamond  domains. 


Simulation  Setup  and  Methodology:  BMI  and  Carbon  Fiber  Interface 

To  calculate  the  interface  thermal  conductance  (ITC)  across  fiber-matrix 
interfaces,  three  combined  systems  were  created  using  the  relaxed  BMI  matrix  and 
different  equilibrated  amorphous  carbon  structures  to  mimic  the  fiber-matrix  interface. 
A  home-grown  code  was  written  in  C++  in  order  to  generate  initial  combined 
geometries  (see  Figure  5(a)).  Within  this  code,  the  BMI  matrix  was  replicated  along 
the  z-dimension  4  times  to  create  a  slab-like  geometry,  and  was  then  combined  with 
the  equilibrated  carbon  surface  structure  to  generate  initial  coordinates  for  modeling 
the  matrix-fiber  interface.  Interaction  between  the  matrix  and  fiber  was  evaluated 
using  short-range  van  der  Waals  interactions^  Here,  s  and  cr  parameters  for  carbon  of 

Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 


the  fiber-surface  were  employed  to  be  0.086  kcal/mol  and  3.4  A,  respectively,  while 
GAFF  force-field  parameters  were  used  for  the  matrix.  A  cut-off  of  12  A  was  used  for 
van  der  Waals  interactions.  The  van  der  Waals  interaction  parameter  cross-terms 
between  the  matrix  components  and  the  fiber  surface  were  calculated  using  Lorentz- 
Berthelot  mixing  rule  [17]. 


Equilibrated  Geometry 


21  nm 


Figure  5.  Representative  snapshot  of  the  (a)  initial  generated  geometry  and  (b)  final  equilibrated 
geometry.  The  color  scheme  for  this  figure  is  same  as  Figure  3  and  Figure  4. 


Once  the  combined  initial  geometry  was  finalized,  NVT  and  NPT  simulations 
were  ran  in  succession  to  equilibrate  the  temperature  and  pressure  of  the  combined 
system,  respectively.  Similar  to  the  previous  BMI  discussion,  while  the  temperature 
was  equilibrated  within  200  ps,  all  the  pressure  components  reached  equilibration  after 
4  ns  of  simulation.  The  testing  for  equilibration  of  NPT  simulations  was  evaluated  by 
confirming  that  all  3  box  dimensions  were  equilibrated  and  were  not  changing  with 
time.  After  equilibration,  non-equilibrium  MD  (NEMD)  simulations  were  perfonned 
to  calculate  interface  thermal  conductance.  The  details  of  the  conductance  simulations 
will  be  discussed  in  an  appropriate  section  for  the  sake  of  continuity  of  discussion. 


RESULTS  AND  DISCUSSION 

Atomic  Partial  Charge  Calculations  using  R.E.D. 

RESP  charge  calculations  were  perfonned  using  RESP-ESP  charge  Derive 
(R.E.D.)  program  as  briefly  mentioned  above.  We  should  point  out  that  RESP  atomic 
charges  are  often  preferred  when  employing  GAFF  force  field  (AMBER)  for  MD 
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simulations.  First,  Materials  Studio®  Visualizer  was  used  to  build  4  different 
molecular  conformations  (models)  of  each  of  the  two  matrix  components.  Then,  all  the 
models  were  optimized  using  HF/6-3 1 G*  level  with  tight  convergence  criteria  in 
GAUSSIAN09  program,  followed  by  molecular  electrostatic  potential  (MEP) 
computations  at  HF/6-3 1G*  level.  Merz  Singh  Kollman  [18]  method  was  used  during 
the  calculation  of  atomic  point  charges  in  Gaussian09.  The  final  electrostatic  charge 
fitting  was  perfonned  using  resp2.2  code  as  provided  with  R.E.D.  program  [11].  The 
output  charges  of  each  component  were  compared  for  different  starting  conformations 
and  were  averaged  out  to  get  final  charges.  We  should  point  out  that  for  both 
monomers,  negligible  differences  were  found  in  the  final  estimated  charges.  Finally, 
the  charge  output  along  with  GAFF  force-field  (for  bonded  and  van  der  Waals 
interactions)  was  imported  into  LAMMPS  fonnat  using  AMBER’s  xleap  program  and 
amber21ammps  python  script  as  provided  with  LAMMPS  [12]. 


Temperature  (K) 

Figure  6.  Plot  of  BMI  matrix  density  as  a  function  of  temperature,  calculated  using  equilibrium  volume 
simulations  (second  set  of  simulations;  see  text).  The  gray  lines  are  fitted  to  the  low  and  high 

temperature  profiles  to  estimate  Tg. 


Thermo-Physical  Properties  of  BMI  Matrix 

To  evaluate  thenno-physical  properties  and  their  variation  with  temperature,  the 
previously  equilibrated  BMI  matrix  (-9500  atoms)  was  heated  to  560  K  in  NPT 
ensemble  for  200  ps.  Thereafter,  two  sets  of  simulations  were  performed  to  model  the 
temperature  behavior  of  different  thermo-physical  properties.  In  the  first  case,  the 
matrix  was  annealed  to  200  K  at  two  different  rates  (0.01  K/ps  and  0.25  K/ps)  from 
560  K.  This  was  done  to  investigate  how  the  cooling  rate  affects  the  glass  transition 
temperature  of  the  matrix.  For  the  second  set  of  simulations,  the  room  temperature 
matrix  was  ramped  up  or  down  to  a  desired  temperature  (between  200  K  and  560  K 
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with  steps  of  40  K)  over  200  ps,  followed  by  1  ns  of  equilibrium  simulations  in  NPT 
ensemble.  In  these  simulations,  the  last  200  ps  of  data  were  used  to  calculate 
equilibrium  volume/density  at  that  temperature. 


Temperature  (K) 

Figure  7:  (Right)  Plot  of  matrix  volume  as  a  function  of  temperature  for  fast  annealing  simulations  (red 
squares),  slow  annealing  simulations  (green  circles)  and  equilibrium  simulations  at  each  temperature 
(blue  triangles).  The  gray  lines  are  fitted  to  the  low  and  high  temperature  profiles  to  estimate  Tg.  The 
green  and  red  vertical  lines  are  shown  to  highlight  the  differences  in  predicted  Tg  due  to  differences  in 

cooling  rates. 


Figure  6  shows  the  variation  of  density  with  respect  to  the  temperature  as 
calculated  from  the  slowest  annealing  simulation  study  (O.OlK/ps).  Based  on  the 
results,  the  room  temperature  density  was  observed  to  be  -1.17  gm/cm  and  matches 
with  the  one  estimated  using  simulations  perfonned  on  the  initial  BMI  structure 
(Figure  3(b)).  Given  the  reported  densities  of  the  individual  monomers  at  room 
temperature  (1.07  g/cm3  for  0,0’  diallyl  bisphenol  A  [19]  and  1.34  g/cm3  for  4,4’- 
bismaleimidodiphenyl  methane  [20]),  In  addition,  it  is  expected  that  the  overall  system 
density  should  increase  with  curing  (cross-linking)  due  to  volume  shrinkage,  as  seen  in 
other  studies  as  well  [21],  and  will  be  investigated  in  the  future.  Thus,  the  value  of 
1.17  gm/cm3  for  un-crosslinked  mixture  is  in  excellent  agreement  with  the  reported 
value  of  1.23  gm/cm  for  cured  matrix  based  on  BMPM  and  DABPA  [22],  assuming 
the  potential  volume  shrinkage  of  the  matrix  during  curing. 

Figure  7  plots  the  volume  as  a  function  of  temperature  and  identifies  the 
importance  of  cooling  rates  on  the  structural  relaxation  of  the  molecules.  Here,  the 
average  volume  for  each  plotted  data  point  was  calculated  using  histogram  binning 
method  with  AT  of  5  Kelvin  for  annealing  simulations.  The  volume  calculation 
strategy  for  equilibrium  simulations  has  been  mentioned  in  a  previous  sub-section 
(refer  to  second  set  of  simulations).  From  the  figure,  it  is  clear  that  a  higher  cooling 
rate  (0.25K/ps;  red  squares)  results  in  larger  system  volume  at  low  temperature  and 
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leads  to  a  higher  value  of  fitted  glass  transition  temperature  ( Tg  ~  490  K).  At  higher 
rates,  the  system  has  lower  chances  to  respond  to  the  non-equilibrium  surrounding 
changes  (such  as  continuous  cooling)  and  thus  leads  to  frustration  in  molecular 
relaxation  [23].  In  such  cases,  the  volume  further  reduces  at  low  temperatures  due  to 
reduction  in  free  volume  and  physical  aging.  On  the  contrary,  relatively  lower  cooling 
rates  (0.01  K/ps;  green  circles)  lead  to  better  equilibrated  structures  and  matches  well 
with  the  data  from  equilibrium  simulations  (blue  triangles). 

It  should  be  pointed  out  that  the  reported  results  for  volume  variation  with 
temperature  and  the  predicted  glass  transition  temperatures  are  for  studied  BMI- 
monomer  mixtures.  It  is  expected  that  after  curing,  the  cross-linked  structure  will 
exhibit  significantly  slower  response  to  cooling  rate  due  to  chain  connectivity  and 
large  molecular  weight,  and  thus  is  expected  to  exhibit  significantly  higher  glass 
transition  temperatures  and  provide  better  comparison  between  simulations  and 
experimental  results  (Tg~  525  K)  [22]. 


Temperature  (K) 

Figure  8:  Plot  of  coefficient  of  linear  thermal  expansion  (in  ppm/°C)  as  a  function  of  temperature  for  the 
BMI  matrix  components,  as  calculated  from  slow  annealing  simulation  run. 


The  coefficient  of  thermal  expansion  is  one  of  many  important  thermodynamic 
properties  in  composite  materials  which  finds  its  use  in  several  problems  of  practical 
importance,  such  as  lattice  mismatch  issues  at  interfaces  between  layers  (films)  of 
different  materials  (composites,  hybrids,  metals,  etc.).  The  coefficient  of  volume 
thermal  expansion  (CVTE)  is  defined  as  [24] 


1  fdV^ 

(X:  =  —  - 

V{dTJp 


(1) 
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where  V,  T,  and  P,  are  the  volume,  temperature,  and  pressure  of  the  system, 
respectively.  For  isotropic  materials,  the  coefficient  of  linear  thennal  expansion 
(CLTE),  /?,  is  related  to  a  as 


p- f  (2) 

Figure  8  shows  the  CLTE  as  a  function  of  temperature  for  the  studied  BMI-5292 
matrix,  as  calculated  from  volume-temperature  data  from  slower  annealing 
simulations.  Here,  a  finite  difference  method  was  used  to  evaluate  /?  at  each 
temperature  T  using  following  equation. 


Pt  = 


Vr-VT 

3xVr  x  (T'-T ) 


(3) 


where,  T’  >  T  and  T'-T  was  fixed  to  be  5  Kelvin  (similar  to  volume  vs.  temperature 
estimation  as  reported  above). 

First,  the  figure  shows  that  CLTE  is  not  constant  with  temperature  but  seems  to  be 
increasing  linearly  (with  fluctuations)  over  the  range  of  300  K  from  10 — 30  ppm/°C, 
although  one  could  argue  the  validity  of  its  linearity  dependence.  Additional 
independent  runs  are  required  to  decrease  the  statistical  fluctuations  and  verify  the 
linear  behavior.  While  these  reported  values  are  smaller  than  reported  values  for 
EPON  based  resins  [21],  the  estimated  values  are  in  fair  agreement  with  CLTEs  of 
other  bismaleimide  resins  based  systems  (Unpublished  work;  please  visit  cited  URLs 
for  comparison  purposes)  [22,  25]. 

Thermal  Conductivity  of  BMI  Matrix 

The  thennal  conductivity  of  the  BMI  matrix  was  calculated  using  NEMD 
simulations  from  a  slab-like  geometry  based  on  Fourier  law  approach  [26].  First,  an 
elongated  slab  was  created  by  replicating  the  equilibrated  BMI-5292  matrix  6  times 
along  z-direction,  followed  by  further  equilibration  in  NPT  ensemble  for  500  ps.  The 
schematic  as  well  as  the  associated  nomenclature  of  NEMD  simulations  for  the 
equilibrated  elongated  slab  is  shown  in  Figure  9. 

As  seen  from  the  figure,  the  outer  edges  of  the  elongated  slab  ( 1  mn)  were  fixed 
and  the  adjacent  ~  2  mn  bins  were  employed  as  hot  and  cold  thennostats,  respectively. 
Two  methodologies,  based  on  constant  temperature  and  constant  energy  flux,  were 
employed  to  calculate  the  thennal  conductivity  as  well  as  investigate  their  potential 
variation  in  predicted  A  values.  For  constant  temperature  techniques  (i.e.,  keeping  the 
hot  and  cold  regions  at  a  certain  temperature),  several  different  heating/cooling 
mechanisms  were  used,  such  as  Nose-Hoover  thennostat,  Berendsen  thennostat  and 
temperature  (velocity)  rescaling  thennostat  as  implemented  in  LAMMPS  package.  For 
these  techniques,  the  hot  region  was  kept  at  350  K  while  the  cold  region  was  kept  at 
250  K  during  the  course  of  NEMD  simulations.  In  order  to  keep  the  regions  at  their 
specified  temperatures,  energy  was  continuously  added  and  taken  off  from  the  hot  and 
cold  regions,  respectively.  NVE  ensemble  was  employed  for  un-thennostated  regions 
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(see  Figure  9).  For  constant  energy  flux  simulations,  a  flux  of  3  kcaFmoFps  was  added 
to  hot  (subtracted  from  cold)  thermostats. 


Hot 

Thermostat 


Un-thermostated  region 


26  nm 


Cold 

Thermostat 


Figure  9.  Schematic  of  NEMD  simulations  performed  on  BMI-matrix  to  calculate  its  thermal 
conductivity  using  different  constant  temperature  and  constant  energy  flux  schemes  as  discussed  in  text. 


Irrespective  of  the  applied  methodology,  a  temperature  gradient  profile  is 
established  along  the  elongated  slab  direction  in  steady-state.  The  temperature  gradient 
can  then  be  calculated  by  dividing  the  slab  into  pre-defined  number  of  small  slabs  with 
equal  thickness  (4  A  for  our  case)  and  calculating  the  temperature  of  each  bin  using 
the  following  equation. 


T;  = 


3Nikg  k=1 


(4) 


where,  A)  is  number  of  atoms  in  ith  slab;  and  m a  and  vk  correspond  to  atomic  mass  and 
velocity  of  atom  k,  respectively. 

For  all  cases,  initial  simulations  were  run  for  about  1  ns  in  order  to  achieve  steady 
state  for  heat  flow.  Once  the  steady-state  was  achieved,  further  simulations  were  run 
for  1.5  ns  for  data  collection.  The  thermal  conductivity  was  then  calculated  using  the 
following  equation 


/ dz 


(5) 


where,  QIAlSt  is  the  heat  flux  across  the  slab  and  clT/dz  is  the  temperature  gradient  in 
steady-state. 

Figure  10  shows  the  cumulative  heat  as  a  function  of  simulation  time  for  different 
investigated  thennostats.  The  initial  non-linearity  is  suggestive  of  the  temperature 
evolution  within  the  system  of  interest.  Once  steady-state  is  reached  (~1  ns),  the 
cumulative  heat  input  becomes  linearly  proportional  to  the  simulation  time.  From  the 
overlapping  data,  the  figure  clearly  suggests  that  there  was  no  noticeable  difference  (< 
5-8%)  in  heat  flux  because  of  different  applied  thennostats.  The  total  heat  flux  values 
(as  obtained  from  the  slope)  are  also  tabulated  in  Table  I  for  relative  comparison. 
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Figure  10.  Plot  of  cumulative  heat  added  into  hot  thermostat  as  function  of  time  for  different  studied 
thermostats.  The  color  scheme  for  different  simulations  is  shown  in  legends. 


Bin  Position  (nm) 

Figure  11.  Plot  of  steady-state  temperature  profde  along  elongated  slab  direction  within  un-thermostated 
region  (See  Figure  9  for  the  schematic  representation).  The  color  scheme  for  different  simulations  is 
shown  in  legends.  All  colors  but  cyan  are  not  readily  visible  due  to  overlapping  of  the  data  points  as  well 

due  to  the  linear  temperature  profde  trend. 
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TABLE  I:  COMPARISON  OF  THERMAL  CONDUCTIVITY  AS  CALCULATED  FROM 
DIFFERENT  INVESTIGATED  METHODS 


Thermostat  Scheme 

Heat  Flux 

(keal/mol/ps) 

Area 

(A2) 

dT/dz 

(K/A) 

Thermal  Conductivity  (A) 
(W/m-K) 

Berendsen 

2.88 

2359.4 

0.46 

0.187±0.002 

Nose-Hoover 

3.12 

2359.4 

0.46 

0.206±0.002 

Temp.  Rescaling 

2.93 

2359.4 

0.44 

0.195±0.002 

Constant  Energy 

3.00 

2359.4 

0.43 

0.203±0.002 

Figure  1 1  plots  the  steady-state  temperature  profile  of  the  BMI-5292  matrix  along 
the  heat  flux  direction.  In  addition  to  the  different  investigated  thermostats,  a 
temperature  profile  from  constant  energy  NEMD  simulations  is  also  plotted.  Similar  to 
Figure  10,  the  estimated  temperature  gradients  (see  Table  I)  were  found  to  be 
indifferent  to  applied  thennostats  (under  same  boundary  conditions  of  250  K  and  350 
K).  Table  I  also  shows  the  estimated  values  of  thermal  conductivity  (using  eq.  2)  for 
BMI-matrix  monomer  systems.  All  the  estimated  values  are  within  -7%  of  each  other. 
The  indifferences  in  estimated  values  due  to  thermostats  are  reflective  of  the  fact  that 
thermostats  at  the  boundaries  do  not  play  a  significant  role  in  determining  thermal 
conductivity  of  the  system.  Such  observations  have  been  reported  before  as  well  for 
silicon  nanowires,  where  a  detailed  investigation  was  carried  out  to  specifically 
observe  these  differences  [27].  It  should  be  noted  that  the  reported  values  are  in  good 
agreement  with  simulations  performed  on  other  ‘amorphous’  epoxy  monomers 
mixtures;  EPON-862  and  DETDA  [28], 

While  it  is  not  pertinent  to  compare  the  estimated  X  values  for  monomer  mixture 
with  X  of  the  ‘cross-linked’  BMI-5292  matrix  [29],  it  is  imperative  to  put  the  estimated 
values  in  perspective.  Previously,  we  have  estimated  that  -90%  crosslinking  of 
amorphous  and  disordered  epoxy  matrix  based  on  EPON-862  and  DETDA  increased 
the  thennal  conductivity  by  -50%  from  their  amorphous  un-cross-linked  counterparts 
[28].  Recently,  it  has  been  reported  that  thermal  conductivity  of  amorphous  polymers 
(polythiophenes)  can  be  significantly  tuned  (up  to  -4.4  W/m-K)  by  increasing  the 
orientational  order,  even  without  order  crystallinity  [30].  Given  the  planar  nature  of 
both  the  monomers  and  the  observation  that  one  of  the  monomers  (BMPM)  exhibit 
crystallinity  (melting  point  156  °C)  [31],  a  similar  increase  could  potentially  be 
possible  for  the  matrix  of  interest.  The  modeling  of  crosslinking  and  the  effect  of 
orientational  order  on  thermal  conductivity  will  be  undertaken  in  the  future. 


Hot  Cold 

Thermostat  Thermostat 

Figure  12:  Schematic  of  NEMD  simulations  performed  to  model  and  predict  interface  thermal 
conductance  at  matrix-fiber  interface  as  discussed  in  text. 
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Interface  Thermal  Resistance  across  Carbon-Fiber/Matrix  Interface 

Figure  12  shows  the  representative  schematic  for  NEMD  simulations  employed 
for  the  calculation  of  interface  thermal  conductance  across  the  BMI-matrix  fiber 
interface.  Three  different  amorphous  forms  (as  retrieved  from  equilibrium  simulations 
of  amorphous  carbon)  of  carbon  fiber  were  used  to  model  the  matrix-fiber  interface. 
We  will  refer  to  them  as  EO,  E2  and  E2,  corresponding  to  cooling  rates  of  1E0  K/ps, 
1E1  K/ps  and  1E2  K/ps,  respectively.  Similar  to  previous  discussion  and  as  seen  in 
Figure  12,  the  outside  edges  were  fixed  for  these  calculations  as  well.  The  adjacent  1.5 
mn  near  both  edges  were  used  as  hot  and  cold  thermostat  regions.  Similar  to  the  BMI 
thermal  conductivity  simulations,  three  difference  thermostats  (Nose-Hoover, 
Berendsen,  and  temperature  (velocity)  rescaling)  were  used  in  hot  and  cold  regions  to 
control  the  temperatures  of  350K  and  250K,  respectively.  As  discussed  prior,  when 
the  energy  is  dumped  in  (or  taken  out)  into  hot  (from  cold)  thermostat,  a  steady-state 
temperature  profile  is  generated  over  time  within  the  system  of  interest.  In  addition, 
near  the  interface,  a  temperature  discontinuity  develops  due  to  mismatch  in  vibrational 
properties  of  the  two  media.  For  our  system  of  interest,  steady  state  was  reached  -  Ins 
after  which  production  runs  of  2  ns  were  run  to  estimate  the  temperature  profile  as 
well  as  temperature  discontinuity  at  the  interface  using  the  binning  method  discussed 
above. 


Bin  Position  (nm) 

Figure  13:  Plot  of  steady-state  temperature  profiles  for  fiber-matrix  composite  systems  (corresponding  to 
different  cooling  rates),  as  employed  for  interface  thermal  conductance  calculations.  The  plot  shows  the 
cases  for  Berendsen  thermostat  simulation  cases.  Other  cases  are  not  shown  for  the  sake  of  clarity.  The 
inset  highlights  the  temperature  decay  in  amorphous  carbon  region  for  better  clarity. 
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The  interface  thennal  conductance  was  then  calculated  using  following  equation. 


A  = 


Q 

AAT 


(6) 


where  A,  A,  and  AT  correspond  to  interface  thennal  conductance,  cross-sectional  area, 
and  temperature  discontinuity  at  the  matrix  fiber-interface.  In  these  simulations,  a 
timestep  of  0.2  fs  was  employed  to  keep  the  total  energy  of  the  system  conserved 
during  the  course  of  the  simulation. 

Figure  13  shows  the  steady-state  temperature  gradient  profiles  for  three  of  the 
nine  studied  systems.  It  can  be  observed  clearly  that:  a)  temperature  gradients  are 
linear  in  matrix  region;  b)  they  are  relatively  less  within  the  carbonaceous  region;  and 
c)  there  exists  a  clear  temperature  discontinuity  at  the  interface  of  the  matrix  and 
carbonaceous  region,  as  discussed  previously. 

The  large  temperature  gradient  within  the  matrix  is  indicative  of  its  low  thennal 
conductivity  and  matches  with  the  calculations  perfonned  on  the  pristine  matrix  as 
discussed  before  in  detail  (please  refer  to  Table  I).  On  the  other  hand,  the  relatively 
lower  slope  (for  the  same  amount  of  heat  flux;  please  see  Table  II  for  heat  flux  values) 
as  seen  from  inset  of  Figure  13  suggests  a  relatively  larger  thennal  conductivity  of  the 
amorphous  matrix.  Based  on  a  linear  fit  to  the  amorphous  carbon  data,  the  thennal 
conductivity  within  the  carbonaceous  region  can  be  estimated  to  be  ~3-4  W/m-K.  We 
should  point  out  that  the  value  range  is  rather  qualitative  due  to  less  number  of  data 
points  in  the  linear  fitting  within  the  amorphous  carbon  region. 


TABLE  II:  OBSERVED  HEAT  FLUX  AND  ASSOCIATED  TEMPERATURE  DROP  VALUES  AS 
STUDIED  FOR  DIFFERENT  COMBINED  FIBER-MATRIX  SYSTEMS  WITH  DIFFERENT 

THERMOSTATS 


Thermostat  Scheme 

Heat  Flux 

(kcal/moFps) 

Area 

(A2) 

AT 

(K) 

Cooling  Rate  1  K/ps 

Nose-Hoover 

3.38 

2728.921 

29.19 

Berendsen 

3.22 

2728.921 

28.22 

Velocity  Rescaling 

3.22 

2728.921 

29.95 

Cooling  Rate  10  K/ps 

Nose-Hoover 

3.28 

2728.921 

30.55 

Berendsen 

3.15 

2728.921 

30.40 

Velocity  Rescaling 

3.46 

2728.921 

31.71 

Cooling  Rate  100  K/ps 

Nose-Hoover 

3.39 

2728.921 

29.65 

Berendsen 

3.31 

2728.921 

29.07 

Velocity  Rescaling 

3.31 

2728.921 

31.01 
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Based  on  the  data  in  Figure  13,  the  linear  region  was  fit  to  calculate  the 
temperature  drop  at  the  interface.  The  temperature  gradients  were  estimated  to  be  ~30 
±  2  K  for  all  studied  cases  (irrespective  of  applied  thennostat  and  nature  of 
carbonaceous  region).  These  values  are  reported  in  Table  II  along  with  the  estimated 
heat  flux  and  cross-sectional  area  of  the  interface.  The  similar  values  of  temperature 
drop  are  indicative  of  the  observation  that  thermostats  do  not  play  an  important  role 
towards  determination  of  thermal  energy  exchange  at  the  interface.  Furthermore,  it 
also  suggests  that  the  energy  exchange  is  predominantly  governed  by  dispersive  van 
der  Waals  interactions.  As  the  local  atomic  structure  of  the  amorphous/nano¬ 
crystalline  carbonaceous  region  is  not  of  noticeable  significance  towards  thermal 
energy  exchange,  as  it  is  the  interface  which  serves  as  a  bottleneck  for  efficient 
thermal  conduction. 


O 

4-» 

£  Thermostat  Scheme 


Figure  14.  Predicted  values  of  interface  thermal  conductance  in  units  ofMW/m2-K  for  all  studied 
combined  matrix-fiber  systems  with  respect  to  different  applied  thermostats.  Color  scheme  is  shown  in 

legends  within  the  figure. 


The  calculated  interface  thermal  conductance  values  (as  calculated  from  eq.  6)  are 
shown  in  Figure  14  for  different  studied  cases.  As  seen  from  the  figure,  for  all  cases, 
the  interface  conductance  values  were  found  be  in  a  narrow  range  of  -26.5 — 29.5 
MW/nT-K  with  a  relative  average  error  of  -4  MW/m'-K.  The  primary  source  of  error 
in  the  investigated  values  arose  from  the  error  in  heat  flux  calculations  and  is  expected 
to  be  reduced  by  running  longer  simulations  [32].  We  should  note  that  the  observed 
values  are  in  good  agreement  with  other  reported  studies  involving  carbon  (CNTs)- 
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matrix  (epoxy,  octane,  polymers)  interfaces  which  interact  through  van  der  Waals 
interactions  [33-38]. 


SUMMARY  AND  CONCLUSIONS 

Using  molecular  dynamics  simulations,  we  have  investigated  the  thenno-physical 
properties  of  BMI-5292  matrix  monomer  mixtures  and  its  interface  properties  with 
amorphous  carbon  fiber.  To  the  best  of  authors’  knowledge,  this  is  the  first 
investigation  to  focus  on  BMI-5292  matrix  system  using  full  atomistic  resolution  and 
is  modeled  using  general  amber  force  field.  To  validate  the  force-field,  several  thenno- 
physical  property  values  were  estimated  and  were  found  to  be  in  good  agreement  with 
the  available  experimental  literature. 

Future  research  will  build  on  this  study  and  will  model  and  investigate 
crosslinking  phenomenon  in  this  matrix  system  and  how  different  thennal  and 
mechanical  properties  are  affected  by  the  degree  of  crosslinking.  In  parallel, 
graphitized  carbon  interfaces  will  be  created  to  model  the  effect  of  degree  of 
graphitization  and  its  orientation  on  the  interface  thennal  conductance  with  cross- 
linked  BMI-5292  matrix. 
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Appendix  B.  Gold  Nanorod  Temperature  Sensors 
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ABSTRACT:  A  technique  is  reported  for  measuring  and  mapping  the  maximum  internal 
temperature  of  a  structural  epoxy  resin  with  high  spatial  resolution  via  the  optically  detected 
shape  transformation  of  embedded  gold  nanorods  (AuNRs).  Spatially  resolved  absorption  spectra 
of  the  nanocomposites  are  used  to  determine  the  frequencies  of  surface  plasmon  resonances. 

From  these  frequencies  the  AuNR  aspect  ratio  is  calculated  using  a  new  analytical  approximation 
for  the  Mie-Gans  scattering  theory,  which  takes  into  account  coincident  changes  in  the  local 
dielectric.  Despite  changes  in  the  chemical  environment,  the  calculated  aspect  ratio  of  the 
embedded  nanorods  is  found  to  decrease  over  time  to  a  steady-state  value  that  depends  linearly 
on  the  temperature  over  the  range  of  100—200  °C.  Thus,  the  optical  absorption  can  be  used  to 
determine  the  maximum  temperature  experienced  at  a  particular  location  when  exposure  times  exceed  the  temperature- 
dependent  relaxation  time.  The  usefulness  of  this  approach  is  demonstrated  by  mapping  the  temperature  of  an  internally  heated 
structural  epoxy  resin  with  10  pm  lateral  spatial  resolution. 

KEYWORDS:  gold  nanorods ,  polymer  nanocomposites,  temperature  sensing  thermal  reshaping  structural  composites 


■  INTRODUCTION 

Gold  nanorods  (AuNRs)  are  known  to  change  their  shape  from 
cylindrical  to  ovoid  or  spheroid  under  the  influence  of  both 
environmental  heat  and  direct  laser  irradiation.  -"  This  thermal 
shape  transformation  leads  to  a  correlation  between  the 
temperature  of  a  gold  nanorod  and  its  optical  properties, 
which  are  dominated  by  shape-dependent  plasmon  resonances. 
The  shape  change  is  primarily  driven  by  the  competition 
between  surface  energy  minimization  across  both  sides  of  the 
interface  and  recrystallization  within  a  surface  melt  layer  (or 
virtual  melt).6'8-11  The  exact  mechanism  of  the  shape 
transformation  of  AuNRs  is  not  fully  understood,  and  no 
analytical  expression  exists  to  describe  the  characteristic 
relaxation  time. 

AuNRs  are  reshaped  to  an  extent  which  depends  on  a 
number  of  factors,  including  the  absolute  size  of  the 
particle,  -  crystal  structure  and  defects,2'11,16’1  and  the 
nature  of  the  interface  with  the  surrounding  medium.  "  ’ 
Because  thermal  reshaping  is  strongly  dependent  on  the 
nanorod’s  environment,  the  rate  and  amplitude  of  the  shape 
change  at  a  particular  temperature  can  vary  dramatically.  For 
example,  AuNRs  that  have  been  stabilized  by  enclosures  of 
metals  or  oxides  can  have  surface  melting  onsets  that  differ  by 
hundreds  of  degrees  from  uncoated  rods.21-  Moreover,  when 
the  surrounding  matrix  properties  change  upon  heating,  the 
correlation  between  temperature  and  spectra  becomes  more 
complicated.  Recently  Liu  et  al.  observed  that,  when  held  at 


elevated  temperatures  for  extended  times,  AuNRs  in  poly- 
(methyl  methacrylate)  (PMMA)  evolve  asymptotically  toward 
a  final  shape  that  is  between  the  as-grown  cylinder  and  the 
thermodynamically  favorable  spheroid.4  Similar  observations 
have  been  made  by  Ng  and  Cheng8  as  well  as  Petrova  et  al.7 
Even  after  many  days  of  heating,  the  rods  in  PMMA 
nanocomposites  did  not  transform  completely  to  spheres.  Liu 
et  al.  demonstrated  the  usefulness  of  the  temperature- 
dependence  of  the  final  nanorod  shape  by  using  a  thermal 
gradient  to  produce  a  color  gradient  in  a  composite  film,4  but 
they  did  not  attempt  to  quantify  the  temperature  gradient  in 
the  film  via  the  nanorod  spectra.  In  this  work,  we  have 
developed  a  technique  that  enables  the  microscale,  ex  situ 
measurement  of  temperature  within  a  polymer  nanocomposite 
by  exploiting  this  thermally  induced  shape  change.  We  calibrate 
the  temperature  response  of  the  rods  by  measuring  their  shape 
after  long  isothermal  holds  while  accounting  for  changes  in  the 
local  dielectric  constant.  We  then  use  this  calibration  to 
determine  the  maximum  temperature  experienced  within  a 
nanocomposite  with  a  lateral  spatial  resolution  of  10  //m.  We 
demonstrate  the  usefulness  of  this  technique  by  applying  it  to 
measure  the  thermal  gradients  inside  a  structural  epoxy  resin 
containing  a  single  6  pm  carbon  fiber  that  is  electrically  heated. 
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Figure  1.  Examples  of  the  spectral  characteristics  necessitating  the  new  analytical  approach  presented  here.  (Left)  SPRpeak  positions  as  a  function  of 
time  for  two  different  isothermal  holds.  The  vertical  dashed  line  represents  the  end  of  the  initial  heating  ramp  period.  The  inset  shows  a  subset  of  the 
absorption  spectra  at  180  °C.  (Right)  Aspect  ratios  as  a  function  of  time  determined  simultaneously  via  SAXS  and  spectroscopic  analysis  for  the 
same  samples. 


The  same  approach  is  extendable  to  much  finer  spatial 
resolutions  via  improved  optical  methods,  and  it  may  be 
expanded  to  larger  temperature  ranges  and  alternative  material 
systems  by  using  different  nanoparticles. 


■  SPECTROSCOPIC  ANALYSIS 


At  ~40  nm  in  length,  the  rods  in  our  AuNR-Epoxy  are  small 
enough  to  treat  with  the  dipole  approximation  (ignoring  the 
size-dependent  dielectric  correction)  without  consideration  of 
quantum  scattering  effects."4-"”  Under  these  conditions,  in  a 
nonabsorbing  medium,  one  can  calculate  the  extinction  cross 
section  spectrum  from  the  aspect  ratio  according  to  Gans’ 
extension  of  the  classical  Mie  theory.  Following  the  treatment 
of  Eustis  and  El-Sayed1N  we  write  the  extinction  y  as  a  function 
of  photon  energy  E  as 
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where  N  is  the  number  of  scatterers,  V  is  the  volume  of  each 
scatterer,  and  el  and  e2  are  the  real  and  imaginary  dielectric 
constants  of  the  scatterer,  respectively.  The  nonabsorbing 
surrounding  medium  is  considered  to  have  a  frequency- 
independent,  real  dielectric  constant  em.  The  Pj  are  three 
depolarization  factors  which  depend  only  on  the  aspect  ratio  a. 
These  factors  lead  to  two  peaks  in  the  absorption  spectrum 
corresponding  to  the  surface  plasmon  resonances  (SPRs) 
parallel  ( EL )  and  perpendicular  ( ET )  to  the  rod  axis.  Others 
have  analyzed  AuNR  absorption  spectra  using  linear  approx¬ 
imations  of  EL(a;  em ),  and  have  ignored  the  small  but 
measurable  changes  in  ET.1S’~  However,  in  complex  systems 
the  changes  in  local  dielectric  can  obfuscate  the  relationship 
between  the  absorption  spectrum  and  the  temperature.  In  an 
epoxy  nanocomposite  (AuNR-Epoxy),  for  example,  the 
dielectric  constant  of  the  matrix  can  change  as  a  result  of 
chemical  reactions  (such  as  cross-linking),  glass  transitions,  and 


photothermal  degradation.  As  demonstrated  in  this  work,  this 
makes  a  direct,  unambiguous  correlation  between  absorption 
peak  position  and  temperature  impossible.  The  collective 
influence  of  a  simultaneously  changing  dielectric  and  aspect 
ratio  is  apparent  in  Figure  1,  which  shows  the  evolution  of  the 
transverse  (TSPR)  and  longitudinal  (LSPR)  resonance  energies 
in  AuNR  Epoxy  for  isothermal  holds  at  140  and  180  °C  (as 
measured  by  the  total  extinction  in  optical  transmission). 

The  SPR  peak  positions  in  Figure  1  are  identified  by  fitting 
each  spectral  region  to  a  Gaussian  line  shape  on  top  of  a  locally 
linear  background.  At  lower  temperatures  the  LSPR  shifts  to 
higher  energies  as  expected,  but  at  temperatures  above  about 
150  °C  the  LSPR  first  blue  shifts  and  then  slightly  red  shifts.  If 
the  longitudinal  peak  were  used  as  a  direct  measure  of  nanorod 
morphology  in  this  sample  then  we  would  conclude  that  the 
nanorods  were  first  contracting  and  then  lengthening.  This  is 
not  consistent  with  the  currently  accepted  physical  picture  of 
the  shape  transformation  process.  To  illustrate  this  departure 
from  expected  behavior  we  have  included  a  dashed  red  line  in 
the  180  °C  curve  corresponding  to  a  simple  exponential  fit  to 
the  data  between  150  and  400  s.  This  line  represents  the  LSPR 
evolution  that  would  be  expected  based  on  the  current 
literature  on  the  subject. 

To  understand  why  the  LSPR  deviates  from  this  expected 
behavior,  it  is  helpful  to  analyze  the  concurrent  changes  in  the 
TSPR  The  TSPR  peak  is  much  broader  than  the  LSPR  peak, 
primarily  because  of  the  presence  of  spherical  particles  in  the 
starting  material,  inhomogeneous  size  distribution,  and 
homogeneous  broadening  from  surface  scattering.  Even  so,  it 
is  possible  to  track  the  very  small  changes  in  the  TSPR  peak 
position  because  the  uncertainty  in  the  fit  parameter  is  small. 
The  TSPR  resonance  in  each  case  increases  sharply  once  the 
target  temperature  is  reached  and  then  begins  to  gradually 
decrease.  While  the  spectral  shifts  are  small  compared  to  the 
peak  width,  the  movement  is  more  apparent  when  the  spectra 
are  stacked  on  top  of  each  other  as  in  the  inset  in  Figure  1. 
Again,  the  nonmonotonic  changes  in  the  TSPR  peak  energy  are 
inconsistent  with  a  model  that  only  considers  changes  in  the 
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nanorod  aspect  ratio.  To  understand  this  phenomenon,  we 
invoke  some  general  properties  of  Mie— Gans  scattering  in  gold 
nanorods  as  described  below.  Referring  to  eq  1,  for  a  particular 
dielectric  environment  the  longitudinal  mode  energy  EL  is  red- 
shifted  and  the  transverse  mode  energy  ET  is  blue-shifted  as  the 
aspect  ratio  increases.  On  the  other  hand,  as  em  increases  both 
El  and  Et  shift  to  lower  energies  (see  Figure  SI  for  an 
illustration  of  this  effect).  The  observed  evolution  of  the  SPR 
energies  reflects  a  material  system  in  which  the  average  aspect 
ratio  of  the  nanorods  monotonically  decreases  but  the  effective 
surrounding  dielectric  does  not  change  monotonically, 
presumably  due  to  chemical  processes  in  the  resin.  This 
behavior  is  consistent  with  the  well-known  cure  temperature 
dependency  of  the  index  of  refraction  in  epoxies.24  At  short 
times,  the  dielectric  constant  of  the  epoxy  decreases  because  of 
thermal  expansion.  At  longer  times,  the  dielectric  constant 
increases  as  the  degree  of  cross-linking  increases,  resulting  in  a 
net  redshift  in  both  SPR  energies. 

To  develop  a  more  robust  mapping  between  the  absorption 
spectrum  and  the  temperature,  we  begin  with  the  ansatz  that 
the  two  physical  quantities  a  and  em  are  uniquely  determined 
by  the  two  SPR  absorption  peaks.  Specifically,  we  use  the 
theoretical  /(E)  for  aspect  ratios  between  1  and  5  and  local 
dielectric  between  1  and  10,  along  with  values  for  eq  and  e2  of 
gold  from  Johnson  and  Christy,'  to  find  the  SPR  peak  energies 
El  and  ET.  We  thereby  associate  each  pair  of  (a;  em)  to  the 
corresponding  (El;  Et)  and  then  fit  the  two  resonances  to  a 
smooth  analytical  function  of  these  two  variables.  We 
approximate  these  relationships  as  EL  =  /L(fl)gL(em)  and  -Et  = 
/tMutCGii)  where/and  g  are  polynomials  of  sufficient  order  to 
achieve  the  desired  precision.  This  can  be  rewritten  as 

Ei(a,  O  =  Y,  cijka’ek 

j’k  (2) 

for  each  of  the  two  resonances  (i  =  L,  T).  These  equations  are 
then  solved  numerically  for  a  given  pair  of  measured  (El;  Et) 
to  simultaneously  determine  a  and  em  from  the  absorption 
spectrum.  For  AuNR-Epoxy,  and  over  a  finite  range  of  aspect 
ratios  and  dielectric  constants,  we  used  4th  order  polynomials 
for  /  and  g  (see  Table  SI  for  function  parameters),  and  the 
calculated  aspect  ratios  agree  well  with  independently  measured 
aspect  ratios  as  discussed  below. 

The  average  aspect  ratios  of  the  nanorods  in  a  subset  of  the 
postheated  AuNR-Epoxy  samples  were  measured  by  determin¬ 
ing  the  length  and  diameter  of  individual  rods  in  TEM 
micrographs.  A  few  representative  images  are  shown  in  the 
Supporting  Information.  Unfortunately,  the  sample  geometry 
and  composition  makes  such  measurements  difficult,  and  the 
dilute  concentration  makes  it  impossible  to  survey  large 
numbers  of  rods  efficiently.  Therefore,  we  used  small-angle 
X-ray  scattering  (SAXS)  to  measure  the  average  aspect  ratio  of 
a  relatively  large  sample  volume.  SAXS  data  were  collected  at 
the  Advanced  Light  Source  both  during  and  after  heating 
various  samples.  Absorption  spectra  were  obtained  simulta¬ 
neously  and  analyzed  independently.  The  agreement  between 
the  two  independent  measurement  techniques,  shown  for  two 
different  temperatures  in  Figure  1,  is  remarkable.  It  provides 
strong  validation  for  the  spectral  analysis  described  in  this  work. 

■  TEMPERATURE  CALIBRATION 

The  AuNR  shape  transformation  is  an  irreversible  thermody¬ 
namic  process  with  an  equilibrium  state  that  is  determined  by 


the  temperature  and  by  surface  and  bulk  energies.  The  rate  of 
change  of  the  aspect  ratio  in  AuNR-Epoxy  is  several  orders  of 
magnitude  smaller  than  the  thermal  motion  of  the  atoms,  and 
at  least  1  order  of  magnitude  smaller  than  the  overall  heating 
rate  of  the  samples.  Under  such  quasi-isothermal  conditions  we 
expect  that  the  final  aspect  ratio  of  the  AuNRs  will  not  be 
strongly  dependent  on  temperature  history.  Thus,  the  final 
aspect  ratio  of  the  AuNRs  after  heating  in  a  given  environment 
is  determined  uniquely  by  the  maximum  steady  state 
temperature  experienced  (see  Figure  S2).  Because  of  changes 
in  the  local  dielectric,  the  energy  of  either  SPR  resonance  alone 
is  not  enough  to  determine  the  temperature.  However,  if  the 
AuNRs  are  exposed  to  elevated  temperatures  for  long  periods, 
then  their  spectroscopically  determined  aspect  ratios  can  be 
used  as  a  measure  of  the  steady  state  temperature  they 
experienced.  We  calibrated  the  AuNR  temperature  response  by 
holding  a  set  of  samples  at  various  uniform  temperatures  and 
recording  the  absorption  spectra  over  time.  Inhomogeneous 
broadening  from  the  size  and  aspect  ratio  distribution  of  our 
AuNRs  leads  to  resonance  peaks  that  are  adequately  fit  by  two 
Gaussian  curves,  for  EL  and  £x,  respectively,  along  with  a 
second-order  polynomial  background.  A  few  of  the  spectra  for  a 
sample  held  at  160  °C  are  shown  in  Figure  2a  along  with  their 
spectral  fits  and  the  time  evolution  of  the  peak  positions 
extracted  from  the  fits  (inset).  Unlike  the  previously  reported 
behavior  of  AuNRs  in  solution,4’  neither  EL(t)  nor  ET(t)  in 
AuNR-Epoxy  are  well  fit  by  any  simple  monotonic  function. 

Once  the  time-dependent  aspect  ratio  aT (t)  is  extracted  from 
each  set  of  spectra  (by  solving  eq  2)  the  relevant  dynamics 
become  more  apparent.  Figure  2b  highlights  these  dynamics  by 
showing  aT(t )  for  T  =  160  °C.  The  time  evolution  is  fit  well  by 
the  biexponential  decay 

aT(t)  =  a1e~t/Zl  +  a2e~t,H  +  ax  (3) 

aoo  =  ao  -  («i  +  a2 )  (4) 

where  a0  is  the  original  aspect  ratio  of  our  AuNRs  and  the  other 
terms  are  fit  parameters  related  to  the  amplitude  (av  a2)  and 
rate  (tv  t2)  of  the  shape  transformation  at  the  temperature  T. 
The  shape  evolution  is  not  well  fit  by  a  single  exponential,  as 
shown  in  Figure  2b.  The  inset  of  Figure  2b  shows  the  calculated 
aspect  ratios  for  a  range  of  temperatures  between  100°  (blue) 
and  200  °C  (red),  each  of  which  is  well  fit  by  the  same 
biexponential  function  with  different  (temperature-dependent) 
fit  parameters.  These  fits  have  R 2  values  of  over  0.98  and 
average  residuals  of  less  than  0.05%. 

The  origin  of  the  biexponential  character  of  the  shape 
evolution  is  not  immediately  obvious.  Much  of  the  previous 
work  has  described  the  process  as  diffusion-limited,  and  some 
have  gone  so  far  as  to  calculate  an  activation  energy  for  the 
shape  transformation.1 1  ~~  Liu  et  al.  proposed  that  the  diffusion 
of  the  AuNR  capping  ligands  into  the  polymer  host  could  play  a 
role  in  destabilizing  the  surface,  resulting  in  shape  trans¬ 
formations  at  temperatures  above  the  glass  transition  temper¬ 
ature  of  the  polymer.4  Unfortunately,  this  does  not  explain  why 
the  AuNRs  in  polymers  evolve  toward  a  shape  that  is  not 
spherical,  even  at  temperatures  far  above  the  glass  transition  of 
the  host.  Some  computational  models  have  indicated  that  the 
reshaping  is  not  diffusion  limited,  but  is  governed  in  part  by  the 
barrier  to  nucleation  of  new  facets. 4  Our  results  clearly  indicate 
that  more  than  one  temperature-dependent  phenomenon 
influences  the  shape  transformation  of  gold  nanorods  in 
polymers.  Still,  it  is  unclear  whether  or  not  the  steady-state 
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Figure  2.  Illustration  of  the  process  for  determining  the  final  aspect 
ratio  using  T  =  160  °C  as  example  data,  (a)  Raw  data  (solid  black)  at 
four  different  times  along  with  Gaussian  fits  (dashed  green)  with 
background  subtracted.  Insets  show  peak  positions  as  a  function  of 
time,  (b)  Aspect  ratio  calculated  from  eq  2  as  a  function  of  time.  The 
biexponential  fit  is  the  solid  red  line  through  the  data  and  the  solid 
blue  curve  represents  the  (poor)  single  exponential  fit.  The  horizontal 
line  at  a  =  3.12  shows  the  asymptotic  limit  of  the  fast  process,  while 
the  lower  line  at  a  =  2.79  shows  the  total  (fast  plus  slow)  limit 
corresponding  to  ax.  The  horizontal  line  at  a0  =  3.85  shows  the 
original  aspect  ratio.  The  inset  shows  the  aspect  ratio  evolution  in 
AuNR-Epoxy  for  a  range  of  temperatures,  comprising  the  calibration 
data. 


aspect  ratio  determined  from  the  fits  actually  corresponds  to  an 
equilibrium  morphology  of  the  nanorods.  Previous  observations 
of  nanorod  spectra  during  longer  heating  periods  suggest  that 
the  rate  of  shape  change  reaches  a  limit,4-  and  the  linear 
region  at  long  times  provides  an  upper  bound  on  the 
uncertainty  associated  with  longer  heating  times.  For  example, 
after  600  s  at  180  °C  the  rate  of  change  ( da/dt )  is  only  1  X 
10-10  s-1.  If  the  evolution  asymptotically  approaches  this  linear 
rate  of  change,  rather  than  a  constant  aspect  ratio,  then  it  would 
take  an  additional  2.4  years  of  heating  at  180  °C  to  see  a  1° 
difference  in  the  calculated  temperature.  Because  the 
asymptotic  aspect  ratio  in  AuNR-Epoxy  is  a  function  of 
temperature  only,  the  temperature  calibration  consists  of  the 
relationship  a0C(T),  which  is  approximately  linear  (R2  =  0.95) 
between  100°  and  200  °C.  In  Figure  3,  we  plot  the  calibration 
data  along  with  the  linear  fit.  It  is  worth  noting  that  this  linear 
fit  extrapolates  to  an  asymptotic  aspect  ratio  of  3.73  at  a 
temperature  of  40  °C  (the  resin  processing  temperature), 
which  is  within  3%  of  the  fit  value  a0  =  3.85.  This  fit  value  for  a0 
is  consistent  within  ±0.05  across  all  of  the  calibration  samples, 


which  provides  further  evidence  of  the  robustness  and  fidelity 
of  our  data  analysis  technique.  The  error  bars  of  ±0.05  in  the 
calculated  aspect  ratios  shown  in  Figure  3d  are  estimated  from 
the  uncertainty  in  the  peak  assignments  from  the  spectral  fits. 
This  is  similar  to  the  RMS  error  of  ±0.053  in  the  linear  fit  of 
and  corresponds  to  an  uncertainty  in  the  calculated  temper¬ 
ature  of  6.6  °C.  If  the  LSPR,  rather  than  the  aspect  ratio,  were 
used  to  measure  the  temperature  then  the  calibration  process  in 
AuNR-Epoxy  would  be  impossible.  No  asymptotic  limit  would 
be  found  for  temperatures  above  around  150  °C  (see  Figure  l), 
and  so  the  temperature  calibration  would  be  ill-defined.  At 
these  temperatures,  the  decrease  in  EL  is  accompanied  by  a 
decrease  in  the  ET,  which  allows  decoupling  of  the  dielectric 
changes  from  the  morphological  changes.  Without  this 
correction,  we  would  calculate  a  temperature  that  is  more 
than  20  °C  lower  than  the  actual  temperature.  This  difference 
would  be  temperature-  and  time-dependent,  and  the  method 
would  predict  temperatures  in  poor  agreement  with  the  FLIR 
measurements  and  COMSOL  calculations  as  discussed  below. 

To  verify  the  accuracy  of  this  new  temperature  measurement 
we  prepared  a  rectangular  slab  of  AuNR-Epoxy  with 
dimensions  8.5  X  37.5  X  2.2  mm  and  heated  it  at  one  end 
with  a  hot  stage  while  cooling  at  the  opposite  end  with  a  copper 
coldfinger  immersed  in  water.  We  recorded  the  surface 
temperature  of  the  sample  using  an  infrared  camera.  After 
heating  for  3  h  the  absorption  was  mapped  point  by  point  using 
the  microspectrophotometer.  A  3  X  10  mm  section  of  the  slab 
near  the  hot  end  is  shown  in  Figure  3a,  b. 

A  quantitative  comparison  between  the  two  methods  is 
obtained  by  averaging  the  temperatures  in  each  row  in  the 
spectral  maps.  We  calculated  the  average  linear  temperature 
profiles  from  the  colder  side  to  the  hot  side  using  both  the 
radiometric  and  spectroscopic  methods  and  the  results  are 
shown  in  Figure  3c.  The  temperature  maps  and  average 
temperature  profiles  are  qualitatively  similar,  with  both 
measurements  indicating  hot  and  cool  areas  as  expected.  The 
average  variance  between  the  spectroscopically  determined 
temperatures  and  a  second  order  interpolation  curve  through 
the  FLIR  temperatures  over  the  measured  range  is  7  °C.  At 
higher  temperatures  (>200  °C)  the  deviation  may  be  due  to  the 
extrapolation  of  a^T)  beyond  the  calibration  range.  At  lower 
temperatures  (<100  °C)  the  stochastic  fluctuations  in  the 
calculated  aspect  ratio  are  large  with  respect  to  the  overall 
change  in  aspect  ratio. 

■  APPLICATION  TO  AN  INTERNALLY  HEATED 
SYSTEM 

Finally,  we  demonstrate  the  utility  of  this  novel  technique  by 
measuring  the  interior  temperature  of  an  internally  heated 
specimen  ex  situ.  We  prepared  a  model  system  consisting  of  a 
single  carbon  fiber  surrounded  by  AuNR-Epoxy  in  a  traditional 
dogbone  configuration  as  shown  in  Figure  4.  This  single  fiber 
dogbone  (SFDB)  was  heated  by  delivering  electrical  power 
directly  to  the  fiber.  We  modeled  the  system  with  COMSOL, 
and  Figure  4d  shows  the  computational  results,  from  a 
stationary  finite  element  method  (FEM)  analysis  in  COMSOL, 
which  agree  qualitatively  with  our  observations  via  optical 
radiometry  (Figure  4c). 

Cross-sections  approximately  1.5  mm  thick  at  different 
positions  along  the  fiber  axis  allow  the  measurement  of  thermal 
gradients  inside  the  SFDB  via  the  method  described  above. 
Examples  of  these  cross  sections  are  shown  in  Figure  5.  In  a 
computer-automated,  point-by-point  algorithm  the  aspect  ratios 
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Figure  3.  Temperature  calibration  for  the  nanorod  temperature  sensors,  (a)  Infrared  image  of  thermal  gradient  across  a  slab  of  AuNR-Epoxy.  (b) 
Spectrally  calculated  temperature  map  of  the  postheated  slab,  (c)  Average  temperature  profile  in  the  direction  of  the  gradient  for  each  image. 
Illustrations  of  AuNRs  with  different  aspect  ratios  show  the  physical  scale  of  their  transformation,  (d)  The  final  aspect  ratio  a  (blue  circles)  from  the 
spectral  calculations  along  with  values  measured  by  TEM  (green  circle  crosses)  and  SAXS  (red  squares).  Insets  show  example  raw  data  for  the  TEM 
and  SAXS  measurements.  The  linear  fit  a00(T)  is  shown  as  a  dashed  line. 


were  determined  from  the  absorption  spectra  and  then 
converted  to  temperature.  Figure  5c  shows  the  spectroscopi¬ 
cally  determined  temperature  maps  for  cross  sections  A  and  B 
from  one  such  sample  (see  Figure  4).  In  the  low-resolution 
raster  scans,  the  large-scale  temperature  profile  is  evident,  but 
the  intrinsic  spatial  averaging  blurs  the  temperature  variations 
near  the  fiber  where  the  temperature  gradient  is  the  highest.  At 
high  resolution,  the  steep  thermal  gradients  near  the  fiber  are 
more  apparent,  but  because  the  fiber  is  not  perfectly 
perpendicular  to  the  sample  surface  the  signal  contains 
contributions  from  the  projection  of  the  fiber  throughout  the 


1.5  mm  sample  thickness.  This  causes  an  apparent  stripe  of 
high  temperatures,  and  the  same  phenomenon  is  responsible 
for  the  blurring  in  the  magnified  visible  light  image  in  Figure  5a. 

To  quantitatively  compare  the  spectral  temperature  maps 
with  the  FEM  model,  we  use  the  temperature  profiles  along  a 
straight  line  through  the  fiber  at  cross  sections  A  and  B.  Figure 
5c  shows  the  temperature  profiles  for  the  AuNR  method  along 
with  the  FEM  predictions  for  the  set  of  parameters  that 
minimizes  the  square  deviation  between  the  model  and  the 
data.  The  fitted  model  parameters  of  total  input  power  (P  = 
1.09  W),  convective  heat  transfer  coefficient  (h  =  11  W/(m2 
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Figure  4.  SFDB  after  Joule  heating  of  the  carbon  fiber,  (a)  True  color 
image.  A  solid  blue  line  is  drawn  to  indicate  the  fiber  position,  which  is 
barely  visible  at  this  magnification.  Thick,  vertical  lines  indicate  cross 
sections  A  and  B.  (b)  Optical  micrograph  with  exaggerated  color 
saturation,  (c)  FLIR  image  during  heating,  (d)  Surface  temperatures 
from  COMSOL  solution. 


K)),  and  epoxy  thermal  conductivity  (fc  =  0.201  W/(m  K)) 
result  in  an  average  deviation  of  6.4  °C  between  the  model 
predictions  and  the  spectrally  calculated  temperature  at  each 
point.  Importantly,  the  same  set  of  parameters  is  used  to  fit  the 
temperature  profiles  at  both  cross  sections.  This  increases  our 
confidence  in  the  accuracy  of  the  computational  model  and  the 
fidelity  of  this  new  temperature  measurement  technique.  It  also 
demonstrates  how  one  might  generate  a  three-dimensional 
thermographic  map  by  serially  sectioning  a  sample  and  creating 
raster  scans  at  each  slice.  We  note  that  most  of  the 
temperatures  in  section  A  exceed  the  calibration  range  and 
are  based  on  an  extrapolation  of  a00(T),  but  they  still  fit  the 
COMSOL  model  well. 

■  DISCUSSION 

Many  other  thermometric  techniques  are  limited  in  their  ability 
to  probe  spatial  temperature  distributions  in  three  dimensions 
at  such  small  scales.  Thermocouples  are  available  in  sizes  of  a 
few  microns  or  less,  but  spatial  mapping  with  thermocouples  is 
cumbersome  and  impractical  at  high  resolution.'1'”  Radio- 
metric  thermography  (infrared  imaging)  is  possible  with  spatial 
resolution  better  than  10  ptm,  but  these  methods  are  sensitive 
only  to  the  sample  surface.”’  '4  Some  specialized  techniques  for 
measuring  thermal  properties  of  materials  at  small  scales  have 
been  developed  or  improved  in  the  past  few  years  including 
time-dependent  thermoreflectance  (TDTR),”’ '  third-harmon¬ 
ic  thermoresistance  (3-ffl), 37,38  and  scanning  thermal  micros¬ 
copy  (SThM). 34,39,40  These  techniques  also  impose  significant 
restrictions  on  sample  geometry  and  surface  conditioning,  and 
quantitative  measurements  depend  on  accurate  models  of 
thermal  transport  in  the  experiments.  Our  nanorod-based 
technique  complements  these  other  nanoscale  thermographic 


Position  (mm) 


Position  (mm) 

Figure  5.  SFDB  cross-sections,  (a)  Optical  micrograph  with 
exaggerated  color  saturation,  (b)  Spectral  thermal  map  at  coarse 
steps,  (c)  Temperature  profile  for  two  different  cross-sections.  Dashed 
lines  are  COMSOL  calculations,  and  points  are  spectroscopically 
determined  temperatures;  insets  show  the  spectral  maps  for  slices  A 
and  B  in  Figure  4a. 


methods  by  combining  high  spatial  resolution,  a  wide  range  of 
continuous  temperature  response,  and  the  ability  to  record 
maximum  steady-state  temperature  for  ex  situ  measurement.  It 
is  particularly  well-suited  to  the  study  of  polymer  nano¬ 
composites  where  the  AuNRs  can  be  well-dispersed  in  a 
transparent  matrix  and  the  dynamics  of  the  thermally  induced 
shape  change  are  relatively  slow.  Thus,  it  provides  an  important 
tool  for  studies  of  micro-  and  nanoscale  thermal  transport  in 
such  materials. 

The  temperature  resolution  of  the  nanosensors  depends  on  a 
number  of  factors.  The  intrinsic  heterogeneity  of  the  as- 
prepared,  commercially  obtained  AuNRs  results  in  broad 
plasmon  peaks  due  to  a  distribution  of  sizes.  This  limits  the 
precision  of  the  plasmon  peak  determination  and  also  gives  rise 
to  a  spread  in  the  temperature  response  for  a  given  population 
of  nanosensors.  This  limitation  would  be  minimized  by  using 
high  purity,  uniform  size  distributions  of  nanorods.  The 
accurate  determination  of  the  SPR  energies  is  vital  to  the 
determination  of  the  temperature,  and  even  with  monodisperse 
particles  the  uncertainty  associated  with  the  measurement  of 
the  plasmon  resonance  energies  will  ultimately  limit  the 
precision  of  the  temperature  measurement.  The  physical  basis 
for  the  observed  transformation  dynamics  remains  unknown, 
and  this  necessitates  a  new  calibration  for  each  material  system 
in  which  this  technique  is  used.  For  example,  although  all  of  the 
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temperature  calibration  curves  are  well-fit  by  a  biexponential 
function,  there  is  a  qualitative  difference  between  the  dynamics 
above  and  below  the  glass  transition  temperature  of  the  epoxy 
(around  140  °C).  This  seems  to  indicate  that  the  relevant  rates 
are  determined  in  part  by  the  stiffness  or  reactivity  of  the 
surrounding  matrix.  The  calibration  curves  may  also  differ  as  a 
function  of  nanorod  encapsulation  or  growth  technique.  We 
expect  that  the  analysis  approach  reported  herein  will  enable 
theoretical  and  computational  advances  by  providing  more 
detailed  information  about  the  nanorod  transformations  in 
various  environments. 

■  CONCLUSIONS 

We  have  shown  that  AuNRs  may  be  used  as  embedded, 
continuously  varying,  irreversible  temperature  sensors  over  the 
temperature  range  of  100—200  °C  with  an  accuracy  of  at  least  6 
°C.  For  steady  state  thermal  profiles  of  sufficient  duration  this 
method  enables  the  ex  situ  determination  of  temperature  with 
spatial  resolution  down  to  at  least  10  fii m.  These  temperatures 
agree  well  with  surface  temperature  readings  using  a 
commercial  infrared  camera  as  well  as  with  the  predictions  of 
FEM  simulations.  In  principle,  the  ultimate  spatial  resolution  is 
limited  only  by  the  density  of  nanorods,  and  the  nature  of  the 
SPR  probe.  Thus,  nanoscale  resolution  may  be  achieved  with 
the  appropriate  confocal  or  near-field  approach. 

Because  of  the  dynamic  nature  of  the  shape  transformation, 
the  accuracy  of  this  technique  is  subject  to  the  condition  that 
the  nanorod  aspect  ratios  have  effectively  reached  steady-state. 
A  more  complete  understanding  of  the  shape  transformation 
process  would  enable  the  use  of  these  nanosensors  under 
nonsteady  state  conditions.  For  example,  if  the  rate  of  shape 
transformation  is  limited  only  by  the  solubility  of  the  AuNR 
capping  ligands  in  the  host  polymer,  as  suggested  by  Liu  et  al., 
then  an  analytical  expression  for  the  aspect  ratio  as  a  function  of 
time  may  be  developed  from  well-established  chemical  kinetics. 
With  appropriate  physical  boundary  conditions  for  surface 
temperature  and/or  heating  times  the  full  temperature  history 
inside  a  nanocomposite  could  be  calculated  from  a  series  of 
absorption  measurements. 

■  EXPERIMENTAL  SECTION 

Sample  Preparation.  AuNRs  with  average  dimensions  of  10  X  41 
nm  were  purchased  commercially  (Nanopartz:  Organic  Gold  Nano- 
rodz  E12— 10— 808),  dispersed  in  isopropyl  alcohol,  and  mixed  with  an 
epoxy  resin  (Momentive:  EPON  828  +  Huntsman:  Jeffamine  D230) 
to  form  a  rigid  nanocomposite  herein  referred  to  as  AuNR-Epoxy.  The 
AuNR-Epoxy  contains  0.004%  gold  by  volume  (0.06%  by  weight), 
corresponding  to  a  number  density  of  approximately  1  X  1013  mL-3.  A 
survey  of  TEM  images  confirms  that  the  AuNRs  are  well-dispersed  in 
the  matrix.  In  order  to  avoid  prematurely  transforming  the  AuNRs  we 
limited  their  exposure  to  elevated  temperatures  before  the  experiment. 
AuNR-Epoxy  samples  were  processed  in  small  batches  (<15  mg)  to 
avoid  large  exotherms  and  cured  at  40  °C  for  6  h.  The  resulting 
nanocomposite  is  a  structurally  rigid  material  with  approximately  60% 
of  the  cross-link  density  of  a  fully  cured  system  as  determined  by 
differential  scanning  calorimetry  (DSC). 

For  the  calibration  samples  we  prepared  2  mm  thick  samples  of 
AuNR-Epoxy  and  heated  each  in  a  hot  stage  (Mettler  FP82HT)  to 
temperatures  between  100°  and  200  °C.  The  target  temperature  was 
typically  attained  within  1  min  of  heating  and  was  stable  to  within  1  °C 
thereafter  for  the  duration  of  the  experiment. 

Absorption  Spectroscopy.  Absorption  spectra  were  obtained 
using  a  microspectrophotometer  (CRAIC  MSP  12l)  in  transmission 
mode.  Spatially  resolved  spectra  are  obtained  using  a  combination  of 
high  magnification  and  square  aperture  collection  window.  AuNR 


densities  were  designed  to  yield  a  maximum  optical  density  at  the 
longitudinal  plasmon  resonance  peak  between  0.4  and  1.0  over  a  2  mm 
path  length.  Coarse  raster  maps  were  produced  with  a  lateral  raster 
step  size  of  100  film,  a  15X  objective,  and  a  camera  collection  aperture 
of  40  X  40  film.  Finer  resolution  maps  were  produced  with  a  10  film 
step  size,  a  36x  objective,  and  a  camera  collection  aperture  of  5  X  5 
film. 

TEM.  Transmission  electron  micrographs  (TEM)  were  obtained  of 
100  nm  thick  microtomed  slices  of  the  AuNR-Epoxy  after  heating 
using  a  Phillips  CM  200.  Aspect  ratios  were  estimated  by  measuring 
the  longest  distances  along  the  length  and  width  of  the  hemispherically 
capped  cylinders.  These  values  were  then  averaged  over  sample 
populations  of  10—20  rods  for  each  temperature,  and  the  errors  are 
estimated  as  the  sample  standard  deviation. 

SAXS.  Small-angle  X-ray  scattering  (SAXS)  measurements  were 
performed  at  the  SAXS/WAXS  beamline  7.3.3  of  the  Advanced  Light 
Source  at  Lawrence  Berkeley  National  Laboratory  both  during  and 
after  heating.  The  scattering  intensity  l(Q)  was  fitted  with  the  Irena  ? 
analysis  package  using  standard  models  for  the  scattering  of  rods  in 
dilute  systems,  where  Q  =  (4;rsin  6)/X  is  the  magnitude  of  the 
scattering  vector.  The  form  factor  for  a  rod  with  log-normal 
distribution  was  used  to  fit  the  data  with  the  following  fit  parameters: 
AuNR  volume  fraction,  breadth  of  the  size  distribution  ( fi  »  0.15), 
aspect  ratio,  and  mean  radius.  Raw  intensity  versus  Q  data  was  fitted 
using  previously  reported  methods.41  The  fitted  aspect  ratio  is  noisy 
due  to  the  limited  lower  Q  range,  but  a  is  calculated  from  the  mean 
radius  by  assuming  volume  conservation  of  individual  nanorods.  The 
uncertainty  is  estimated  from  the  fitted  radius  distribution. 

FLIR.  Surface  temperature  measurements  were  obtained  with  a 
FLIR  SC620  using  a  40  mm  lens.  The  camera  was  calibrated  by  the 
manufacturer,  and  an  emissivity  of  0.88  was  assumed  for  the  AuNR- 
Epoxy  surface. 

SFDB.  The  single-fiber  dogbones  (SFDB)  used  in  the  internal 
spatial  temperature  maps  comprise  individual  carbon  fibers  (Hexel 
HexTow  IM7)  encased  in  AuNR-Epoxy.  These  were  heated  with  a 
Keithley  2100  source-meter  via  silver  paint  contacts  directly  to  the 
fiber  at  the  ends  of  the  SFDB. 
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Several  of  the  works  cited  in  the  main  text  have  investigated  the  degree  to  which  expres¬ 
sions  such  as  Equation  1  correlate  with  observed  absorption  spectra  for  various  values  of  a 
and  em,  as  well  as  for  subtle  differences  in  e1  and  e2.  We  make  no  specific  claims  about  the 
accuracy  of  this  expression  other  than  to  highlight  the  agreement  of  our  calculated  aspect 
ratios  with  independently  measured  values  from  TEM  and  SAXS.  However,  our  general  ap¬ 
proach  to  the  temperature  calibration  is  motivated  by  the  fact  that  changes  in  a  and  em 
affect  the  predicted  resonance  energies  in  qualitatively  different  ways.  Figure  S-l  illustrates 
this  principle  by  showing  predicted  extinction  coefficients  for  various  values  of  these  two 
parameters. 
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Figure  S-l:  Theoretical  dependence  of  the  scattering  spectrum  for  AuNRs  on  aspect  ratio 
a  and  local  dielectric  m.  The  upper  plots  show  the  extinction  for  aspect  ratios  of  2.5  (solid 
green),  3.0  (dashed  red),  and  3.5  (dotted  blue)  while  the  local  dielectric  is  fixed  at  3.0.  In  the 
lower  plots  the  dielectric  takes  the  same  sequence  of  values  (2.5,  3.0,  3.5)  while  the  aspect 
ratio  is  fixed  at  3.0.  The  ET  resonances  for  each  plot  are  magnified  in  the  insets.  Arrows 
indicate  the  direction  that  the  peak  shifts  as  the  parameters  increase. 
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The  choice  of  a  separated-variable,  4th  order  polynomial  as  the  numerically  solvable 
approximation  for  the  two  Ei(a,em )  (Equation  2  in  the  text)  is  arbitrary  but  justifiable. 
Table  S-l  shows  the  coefficients  Cjk  that  minimize  R 2  between  the  invertible  polynomial 
function  and  Equation  1,  which  is  non-invertible  due  to  the  empirical  values  for  e±  and 
62-  While  the  zero  and  first  order  terms  obviously  dominate,  the  higher  order  terms  are 
not  negligible.  Only  by  including  up  to  4th  order  terms  is  the  RMS  variance  between  the 
polynomial  function  and  the  full  expression  reduced  to  less  than  one  percent  over  the  sampled 
domain  (a  £  [l,5];em  £  [1,10]). 

Table  S-l:  Coefficients  for  the  polynomial  approximation  of  E(a,em )  in  Equation  4.  Index 
j  is  the  row  and  index  k  is  the  column. 


El 


j\k 

0 

1 

2 

3 

4 

0 

1.37017 

1.13551 

-0.31381 

0.03685 

-0.00161 

1 

2.02695 

-1.82677 

0.45046 

-0.05196 

0.00225 

2 

-1.03890 

0.79840 

-0.20623 

0.02459 

-0.00108 

3 

0.21013 

-0.15923 

0.04281 

-0.00519 

0.00023 

4 

-0.01558 

0.01192 

-0.00328 

0.00040 

-0.00002 

Et 


j\k 

0 

1 

2 

3 

4 

0 

2.61635 

-0.25875 

0.03873 

-0.00447 

0.00020 

1 

-0.26946 

0.36822 

-0.09556 

0.01129 

-0.00048 

2 

0.11918 

-0.15128 

0.04150 

-0.00499 

0.00021 

3 

-0.01997 

0.02520 

-0.00709 

0.00086 

-0.00004 

4 

0.00121 

-0.00154 

0.00044 

-0.00005 

0.000002 

We  also  tried  an  alternate  approach  by  employing  both  polynomial  fits  and  numeric 
interpolation  functions  for  em,  but  the  subsequent  inversion  procedure  was  not  robust  over 
the  set  of  experimental  data  in  the  temperature  calibration  procedure  (ie,  numerically  solving 
Equation  2  proved  unstable). 

In  order  to  test  the  uniqueness  of  the  final  shape  at  a  particular  temperature  we  subjected 
several  samples  to  various  thermal  histories.  Figure  S-2  shows  the  time  evolution  of  the 
calculated  aspect  ratio  for  three  of  these  measurements,  each  of  which  was  initially  heated 
to  120°  C  for  one  hour  and  then  allowed  to  cool  to  room  temperature.  Each  sample  was 
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then  reheated  to  a  steady  state  (~  3  hours)  temperature  of  (a)  100°  C,  (b)  120°  C,  and  (c) 
160°  C.  The  final  aspect  ratio  in  each  of  these  cases  is  uniquely  determined  by  the  maximum 
steady  state  temperature  experienced  regardless  of  the  order  of  heating.  For  example,  sample 
(a)  does  not  change  substantially  during  the  second  heating  cycle  because  it  has  already 
experienced  a  higher  temperature.  Sample  (b)  continues  the  same  biexponential  evolution 
as  if  it  had  not  been  interrupted.  The  evolution  of  sample  (c)  changes  character  completely 
during  the  second  cycle  and  asymptotically  approaches  the  same  value  as  the  calibration 
sample  held  at  160°  C.  We  repeated  these  experiments  for  several  temperature  profiles  and 
in  different  resin  systems  and  found  similar  behavior.  Note  that  the  plots  shown  in  Figure  S- 
2  pertain  to  a  different  batch  of  AuNRs  than  those  of  Figure  2  in  the  text;  therefore  the 
specific  asymptotic  values  are  different. 


Figure  S-2:  Aspect  ratio  evolution  for  three  different  samples  initially  heated  and  held  at 
120°  C  for  30  min  and  subsequently  heated  to  100°,  120°,  and  160°  C  respectively.  In  each 
case  the  final  asymptotic  aspect  ratio  is  limited  by  the  maximum  temperature  experienced 
over  the  course  of  the  thermal  history,  irrespective  of  the  time  between  heating  cycles. 


Figure  S-3  shows  a  few  representative  TEM  micrographs  of  AuNRs  after  heating.  The 
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AuNRs  are  dilutely  dispersed  in  the  matrix,  and  so  only  one  rod  at  a  time  is  visible  at 
high  magnification.  We  measured  approximately  100  individual  nanorods  in  this  manner 
and  report  the  measured  aspect  ratios  for  the  data  in  Figure  2  using  the  variance  in  our 
measurements  for  a  particular  temperature  as  the  uncertainty.  However,  it  is  difficult  to 
know,  without  time-consuming  tilting  measurements,  whether  or  not  the  nanorod  axes  are 
perpendicular  to  the  field  of  view.  This  is  especially  true  because  of  the  intrinsic  inhomogene¬ 
ity  in  the  sample  population.  The  as-received  nanorods  already  had  a  noticeable  variation 
in  radius,  so  that  direct  comparison  of  micrographs  between  temperatures  is  ambiguous. 
Therefore,  we  kept  only  the  aspect  ratios  in  the  top  25%  of  all  measurements  at  a  particular 
temperature,  potentially  biasing  the  results. 


Figure  S-3:  Sample  TEM  images  of  AuNRs  in  Epoxy  matrix  after  heating  for  several  hours 
at  the  temperatures  shown. 
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Appendix  C.  COMSOL®  Multiphysics 
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The  effective  optical  and  thermal  properties  as  well  as  the  opto-thermal  response  of 
composite  materials  are  investigated.  A  stochastic  method  is  developed  which  implements 
a  random  microstructure  generation  (RMG)  algorithm  to  produce  a  synthetic  microstruc¬ 
ture  representing  the  geometry  of  unidirectional  continuous  fiber  reinforced  polymer  matrix 
composites  (PMCs).  Simulation  of  the  synthetic  microstructure  determines  effective  ther¬ 
mal  properties,  effective  optical  properties,  and  the  thermal  response  to  optical  energy. 
The  micromechanics  framework  developed  herein  is  also  applied  to  unit  cell  analysis  and 
compared  with  analytical  solutions.  The  metrics  of  effective  transverse  thermal  conductiv¬ 
ity,  effective  attenuation  coefficient,  and  temperature  response  are  studied  as  a  function  of 
fiber  volume  fraction.  We  show  that  some  composites’  effective  thermal  properties  can  be 
adequately  predicted  using  unit  cell  analysis  or  analytic  techniques;  however,  effective  op¬ 
tical  properties  and  thermal  response  to  optical  input  depend  strongly  the  specific  random 
arrangement  of  the  fibers. 


Nomenclature 


a 

e 

P 

A 

Cp 

H 

I 

lo 

L 
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Attenuation  coefficient,  1/m 
Surface  emissivity 
Density,  kg/m3 
Power  absorption  coefficient 
Specific  heat  capacity,  J/(kq  *  K) 
RVE  Height,  m 
Intensity,  W/m 2 
Incident  intensity,  W/m 2 
RVE  Length,  to 
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Q  src 

Heat  Source,  W/m3 

R 

Power  reflection  coefficient 

T 

Temperature,  K  or  °C 

Tr 

Power  transmission  Coefficient 

h 

Heat  transfer  coefficient,  W /(K  *  m2) 

k 

Thermal  conductivity,  W/(K  *  to) 

<7 

Heat  Flux,  W/m2 

rf 

Fiber  Radius,  to 

t 

Time,  s 

Vf 

Subscripts 

Fiber  volume  fraction 

eff 

Effective  property 

f 

Fiber 

m 

Matrix 

normal 

Direction  normal  to  the  surface 

x,V,z 

x-component,  y-component,  z-component 

I.  Introduction 

Polymer  matrix  composites  (PMCs)  continue  to  increase  as  the  material  of  choice  in  modern  aerospace 
systems  on  account  of  the  high  specific  stiffness,  high  specific  strength,  enhanced  environmental  stability, 
and  the  ability  to  tailor  anisotropic  properties.  The  methods  and  processes  to  implement  PMCs  in  aerospace 
systems  are  fairly  well  established;  however  there  remains  ongoing  development  in  the  methods  to  manufac¬ 
ture  and  finish  PMCs.  High  power  lasers  in  particular  have  been  utilized  to  great  effect  in  various  composite 
processing  and  evaluation  applications,  including  trim  operations,  paint  removal,  thermoplastic  consolida¬ 
tion,  laser  machining,  curing,  medical/dental  applications,  and  non-destructive  evaluation  (NDE).1-'  The 
laser  is  utilized  as  a  focused,  variable  intensity  heat  source  in  each  of  these  applications,  which  requires  an 
understanding  of  how  optical  energy  is  transduced  into  heat  and  how  the  heat  is  subsequently  transferred. 

Transducing  optical  energy  to  thermal  energy  within  composites  must  be  well  understood  in  order  to 
optimize  materials  and  processes  associated  with  these  applications;  thus,  we  study  the  interaction  between 
incident  optical  energy  and  thermal  transport  in  the  microstructure  of  PMCs.  Several  analytic  models  can 
predict  the  effective  properties  of  composites  but  they  are  not  valid  or  well  developed  for  all  composite 
systems  across  all  types  of  energy  input  mechanisms,  especially  those  with  high  fiber  volume  fractions,  fiber 
coatings,  finite  thermal  resistances,  and  for  optical  properties. 

Traditional  composite  mechanics  generally  utilize  representative  volume  elements  with  patterned  pack¬ 
ing  or  unit  cells  to  homogenize  the  microstructure  into  an  effective  composite  material  and  then  apply 
laminated  plate  or  similar  theories  to  understand  a  structure.  Stress  analysis  utilizing  the  homogenized 
microstructure  is  highly  accurate;  however  the  predictions  are  not  nearly  as  accurate  for  estimating  other 
effective  properties  such  as  thermal  conductivity,  electrostatic  coefficients,  or  piezoelectric  coefficients.  The 
micromechanics  community  has  adapted  multi-inclusion  models  to  account  for  discrepancies  when  modeling 
multiphase  composites.  Thus  far  multi-inclusion  models  have  been  applied  to  mechanical,  piezoelectric,  and 
thermal  problems  but  not  optical  to  thermal  transduction. 

Analytic  models  are  able  to  predict  the  effective  thermal  conductivity  of  composite  materials,  including 
those  with  interface  resistance  and  fiber  coatings,  but  are  subject  to  certain  restrictions. 8-12  Many  analytic 
models  rely  on  the  simple  rule  of  mixtures  (SROM),  enhanced  rule  of  mixtures  (EROM),  or  the  equivalent 
inclusion  method  (EIM)  to  predict  effective  properties.  Analytic  models  predicting  thermal  conductivity 
have  shortcomings  which  can  be  overcome  using  microscale  finite  element  analysis  (FEA),  such  as  inaccurate 
conductivity  prediction  at  high  volume  fraction  due  to  neglecting  conductive  paths.13,14 

Sihn  and  Roy12  have  explored  the  applicability  of  the  EROM  and  EIM  to  finding  the  transverse  thermal 
conductivity  of  continuous  fiber  composites;  however,  the  EROM  does  not  include  the  effects  of  coatings  or 
interface  thermal  resistance.  Hasselman  and  Johnson9  have  also  developed  a  closed  form  analytic  solution 
based  on  the  equivalent  inclusion  method  (EIM).  Unit  cell  FEA  can  typically  account  for  high  fiber  volume 
fraction,  fiber  coatings,  and  interface  resistance,  but  still  neglects  conductive  paths  that  may  develop  in  a  real 
microstructure.  Because  of  the  limitations  of  analytic  thermal  models  and  unit  cell  methods  to  account  for 
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these  situations  and  the  lack  of  optical  models  developed  for  fiber  reinforced  composites,  a  micromechanical 
approach  with  RMG  is  implemented.  Though  the  micromechanics  approach  is  capable  of  including  fiber 
coatings  and  interface  resistance,  these  effects  are  not  considered  in  this  work. 

Due  to  the  complex  scattering,  absorption,  and  reflection  occurring  via  inclusion  interaction,  unit  cell 
analysis  and  various  rule  of  mixtures  methods  will  not  be  sufficient  to  determine  effective  optical  properties. 
Analytic  prediction  methods  for  effective  optical  properties  of  continuous  fiber  reinforced  composites  are 
considerably  less  developed  than  the  thermal  conductivity  models,  likely  due  to  their  lack  of  application 
as  optical  components.  Some  simplified  models  can  be  used  to  predict  PMC  optical  behavior  such  as 
the  Kubelka-Munk  two  flux  and  four  flux  models,  but  these  require  experimental  evaluation  to  find  fit 
parameters.15  Geometric  optics  or  ray  tracing  has  also  been  used  successfully  to  study  the  optical  properties 
of  unidirectional  composite  systems  where  the  relative  length  scales  and  wavelengths  fall  within  appropriate 
margins.16  In  this  work,  we  demonstrate  the  implementation  of  a  random  microstructure  generation  (RMG) 
algorithm  and  link  optical  interaction  to  thermal  effects  to  enable  the  accurate  prediction  of  effective  optical 
and  thermal  properties,  as  well  as  the  thermal  response  to  optical  energy,  in  a  stochastic  manner.  The 
approach  is  verified,  through  comparison  to  existing  analytical  models  and  literature  values,  for  predicting 
the  effective  thermal  conductivity  of  composites.  The  full  wave  solution  of  Maxwell’s  equations  is  employed 
to  account  for  diffraction  and  interference  effects  that  become  significant  when  the  wavelength  of  the  incident 
radiation  is  similar  in  size  or  larger  than  the  inclusions  and  other  geometric  scattering  features. 


II.  Method 


A.  Synthetic  Microstructure  Generation 


The  challenges  of  high  volume  fraction,  conductive  paths,  and  unit  cells  are  addressed  through  the  use  of  an 
RMG  algorithm  developed  by  Melro  et  al. 1 '  This  RMG  algorithm  is  further  modified  by  Romanov  et  al.18 
to  include  a  random  minimum  fiber  spacing  to  prevent  an  unnaturally  large  amount  of  touching  fibers  and 
verified  this  modification  through  image  analysis  and  micromechanics  simulations;  thus,  the  RMG  algorithm 
with  Romanov’s  modification  is  implemented  in  Matlab  and  used  in  this  work.  The  RMG  code  uses  a 
random,  hard-core  packing  algorithm  to  place  circular  fibers  in  a  rectangular  matrix  until  either  the  desired 
fiber  volume  fraction  ( Vf )  or  the  jamming  limit  (~  50%)  is  reached.  If  this  limit  is  reached,  a  series  of 
heuristic  re-ordering  steps  moves  the  existing  fibers  closer  together  which  allows  new  fibers  to  be  inserted 
in  the  matrix  rich  regions.  An  example  of  a  microstructure  subject  to  the  jamming  limit  and  one  that 
overcomes  this  limit  is  shown  in  Figure  1. 


Figure  1.  The  microstructure  on  the  left  is  subject  to  the  jamming  limit  of  approximately  50%  while  the 
microstructure  on  the  right  has  achieved  a  volume  fraction  of  65%  through  heuristic  reordering 

The  theoretical  maximum  volume  fraction  «  90.69%)  is  achieved  with  perfect  hexagonal  packing. 
Theoretically,  the  heuristic  reordering  algorithm  will  eventually  achieve  this  volume  fraction  if  the  minimum 
fiber  spacing  condition  is  set  to  zero  and  the  algorithm  is  allowed  to  execute  without  iteration  limits.  However, 
the  practical  upper  limit  on  volume  fraction  is  about  70%  for  both  the  RMG  algorithm  and  actual  composites. 

B.  Thermal  Conductivity 

Effective  thermal  conductivity  of  continuous  fiber  reinforced  PMCs  can  be  found  in  a  number  of  ways.  Using 
the  simple  rule  of  mixtures  (SROM)  as  a  parallel  resistive  network  calculates  the  lower  limit  to  transverse 
thermal  conductivity  (assuming  perfect  interfaces)  and  provides  a  good  approximation  of  the  longitudinal 
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thermal  conductivity.13  The  longitudinal  conductivity  approximation  by  the  SROM  assumes  the  fibers  are 
perfectly  aligned  and  that  fiber  waviness  does  not  significantly  affect  the  effective  thermal  conductivity. 
Other  analytic  methods  have  been  used  to  determine  the  effective  thermal  conductivity,  but  only  the  EIM 
and  SROM  are  included  in  this  work  because  others  have  successfully  used  the  EIM  whereas  SROM  is  the 
simplest  model  available.9, 12  The  boundary  conditions  used  for  predicting  the  effective  thermal  conductivity 
via  FEA  are  shown  in  Figure  2. 


Figure  2.  Boundary  conditions  for  predicting  transverse  thermal  conductivity  in  a  continuous  fiber  reinforced 
composite 


The  boundary  conditions  in  Figure  2  are  chosen  based  on  the  work  of  Pramla  et  al.19  which  explores 
and  evaluates  several  methods  of  applying  boundary  conditions  to  an  RVE  to  approximate  the  thermal 
conductivity  of  composites.  Their  work  identified  the  boundary  conditions  shown  in  Figure  2  as  the  most 
accurate  method  of  predicting  effective  thermal  conductivity.  Once  the  steady  state  temperature  field  is 
solved  for  a  given  microstructure,  the  effective  transverse  thermal  conductivity  is  calculated  using  Equation  1. 


Lf0H  hx(y)dy 
HAT 


Keff  = 


(1) 


Where  Kef  /  is  the  effective  thermal  conductivity,  L  is  the  length  of  the  RVE,  H  is  the  height  of  the  RVE, 
AT  is  the  temperature  difference,  and  hx  is  the  heat  flux  in  the  x  direction.  Unlike  analytic  solutions,  very  few 
composite  property  restrictions,  such  as  high  thermal  conductivity  mismatch  ratio,  exist  in  the  FE  solution. 
Unfortunately,  the  size  of  the  RVE  must  be  at  least  approximately  50x  the  average  fiber  radius  (r^)  in  any 
direction  so  as  to  avoid  boundary  effects  and  find  a  good  approximation  for  the  average  heat  flux.19  The  FE 
simulation  allows  us  to  assess  the  variability  due  to  the  stochastic  nature  of  the  microstructure  because  each 
synthetic  generation  has  a  unique  geometry.  Stochastic  variation  is  generally  small  with  respect  to  effective 
thermal  conductivity  when  perfect  interfaces  are  assumed  (no  coatings  or  thermal  interface  resistance). 


C.  Optical  Properties 

Effective  optical  properties  of  random  microstructure  composites  depend  strongly  on  the  arrangement  of 
fibers  due  to  diffraction  and  interference  effects  that  occur  at  fiber  edges  and  between  fibers.  The  fiber  size 
relative  to  the  wavelength  and  the  inter-fiber  spacing  perpendicular  to  wave  propagation  has  a  significant 
impact  on  how  energy  moves  through  the  composite  due  to  absorption,  slit  diffraction  between  fibers,  and 
edge  diffraction  occurring  at  every  fiber  subject  to  irradiation.  Epoxies  and  other  polymers  used  as  matrix 
material  for  PMCs  are  generally  transparent  in  the  near  infrared  spectrum  (700  nm  to  2500  nm)  while  the 
inclusions  (fibers)  may  have  a  wide  range  of  properties  in  this  spectrum.  Carbon  fibers  are  typically  highly 
absorbing  in  the  near  IR  spectrum  while  glass  fibers  are  mostly  transparent  and  other  fibers  such  as  aramids 
may  have  optical  properties  somewhere  in  between.  Due  to  the  complexity  of  arrangements  and  size  of 
the  fibers  and  fiber  spacing  relative  to  the  wavelength,  the  full  wave  solution  is  solved  to  account  for  all 
the  contributing  effects.  In  cases  where  the  wavelength  is  much  smaller  than  the  smallest  feature  size,  the 
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geometric  optics  approximation  known  as  ray  tracing  may  be  used.  The  geometric  optics  (i.  e.  ray  tracing) 
approximation  may  be  appropriate  when  the  fibers  are  highly  absorbing  and  densely  packed  enough  such 
that  diffraction  effects  are  negligible.  For  lower  volume  fractions,  diffraction  effects  become  more  significant 
due  to  the  prevalence  of  fiber  edges  and  gaps  which  further  motivates  the  full  wave  solution  compared  to  the 
geometric  optics  approach  which  neglects  the  enhanced  absorption  and  scattering  caused  by  diffraction. 

The  effective  optical  properties  of  the  composite  are  quantified  through  a  method  similar  to  the  effective 
thermal  properties:  FE  simulation  of  an  RVE.  Finding  the  effective  optical  properties  of  a  composite  first 
uses  the  RMG  algorithm  to  generate  a  random  microstructure.  A  periodic  port  produces  a  plane  wave  which 
propagates  toward  the  synthetic  microstructure.  The  plane  wave  can  be  either  “s”  or  “p”  polarized  with 
respect  to  the  fiber  surface,  but  in  this  work  the  EM  waves  are  always  “p”  polarized  so  that  the  initial 
electric  field  is  perpendicular  to  the  fiber  axis.  To  capture  and  remove  any  portion  of  the  incident  energy 
that  is  reflected  or  transmitted,  perfectly  matched  layers  (PMLs)  surround  the  microstructure  and  the  port 
at  the  right  and  left  boundaries  as  shown  in  Figure  3.  The  PML  boundaries  quickly  attenuate  any  incident 
wave  and  prevent  reflection  of  waves  in  nearly  any  orientation  without  needing  to  have  prior  knowledge  of 
the  polarization  or  propagation  constant.  The  top  and  bottom  boundaries  in  Figure  3  are  perfect  electric 
conductors  (PEC)  which  causes  perfect  reflection  and  simulates  a  semi-infinite  condition.  A  summary  of 
these  boundary  conditions  is  shown  in  Figure  3. 


Figure  3.  Boundary  conditions  and  model  configuration  for  random  fiber  optical  simulation. 

From  the  optical  simulation  shown  in  Figure  3,  the  effective  values  for  power  absorption  (^4),  power  trans- 
mission  (Tr),  power  reflection  (R),  and  attenuation  coefficient  (ae//)  can  be  found.  Tr  is  found  by  integrating 
the  time  averaged  power  density  on  the  right  side  of  the  RVE  whereas  A  is  calculated  by  integrating  the 
power  dissipation  density  within  the  RVE,  or  just  within  the  fibers  if  the  matrix  is  considered  transparent. 
According  to  conservation  of  energy,  all  of  the  incident  optical  energy  must  either  be  transmitted,  reflected, 
or  absorbed.  Thus,  R  is  calculated  by  subtracting  A  and  Tr  from  unity. 

Due  to  the  finite,  random  nature  of  the  RVE,  the  transmission  may  be  non-zero.  This  is  especially  true  for 
low  volume  fractions  of  highly  absorbing  fibers  any  volume  fraction  of  transparent  or  weakly  absorbing  fibers. 
In  order  to  accurately  predict  effective  optical  properties  for  a  given  microstructure,  the  transmission  term  in 
the  RVE  should  be  as  close  to  zero  as  possible  in  order  to  be  considered  optically  thick.  Achieving  an  optically 
thick  RVE  with  unit  cells  is  difficult  at  low  volume  fractions,  thus  the  unit  cells  are  repeated  along  the  wave 
propagation  direction  to  address  this  issue.  For  highly  absorbing  fibers  at  high  fiber  volume  fractions,  the 
depth  at  which  the  transmission  term  approaches  zero  may  be  within  a  few  fiber  diameters.  However  if 
the  fibers  are  transparent,  weakly  absorbing,  reflective,  or  the  volume  fraction  is  low,  the  composite  could 
conceivably  have  a  non-zero  transmission  term  depending  on  the  total  thickness.  Thus,  the  actual  thickness 
of  the  part  must  be  considered  when  finding  the  effective  optical  properties  of  thin  materials  or  transparent 
or  translucent  materials.  Figure  4  shows  the  various  types  of  unit  cell  configurations  and  an  example  of 
the  random  microstructure  applied  to  the  optical  simulation.  Figure  5  highlights  the  difficulty  of  producing 
an  optically  thick  RVE  via  the  unit  cell  methods  shown  in  Figure  4  by  plotting  the  required  depth  to  be 
considered  optically  thick  as  a  function  of  fiber  volume  fraction. 

The  value  for  ae//  is  found  by  first  reducing  the  two  dimensional  power  density  (i.e.  intensity)  to  a  one 
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Figure  4.  Different  types  of  RVEs  showing  electric  field  propagation  and  location  of  the  input  port  (red  line) 
with  the  wave  always  initially  propagating  normal  to  the  port.  All  RVEs  shown  here  are  at  30%  fiber  volume 
fraction  and  100  |im  (20  fiber  radii,  20 rj)  deep. 


Figure  5.  Optical  simulation  depth  required  to  be  considered  optically  thick  for  various  microstructure  gen¬ 
eration  types.  Values  larger  than  the  part  thickness  are  estimated  based  on  effective  attenuation  coefficient. 


dimensional  spatial  average  so  that  power  is  represented  as  a  function  of  depth.  The  attenuation  coefficient 
(aeff)  is  then  found  by  fitting  the  intensity  as  a  function  of  depth  to  the  Beer-Lambert  law  (Equation  2) 
and  extracting  the  exponential  coefficient.15  As  shown  in  Figure  6,  the  Beer-Lambert  law  describes  the 
power  attenuation  well  in  the  RMG  microstructure  which  allows  ae//  to  be  found  using  a  least  squares 
fit.  Analytically,  the  negative  derivative  of  the  Beer-Lambert  law  can  be  used  to  find  the  heat  source  as  a 
function  of  depth  ( Qsrc )  as  shown  in  Equation  3. 

/( x)  =  I0e-a°ff*x  (2) 

Qsrc(x)  =  =  aeffI °e~aeff*x  (3) 

The  left  plot  in  Figure  6  shows  the  FEA  results  for  the  electric  field  propagation  through  a  30%  volume 
fraction  composite  microstructure  generated  using  the  RMG  method.  The  middle  plot  in  Figure  6  shows 
Equation  2  fit  to  the  intensity  as  a  function  of  depth.  Of  the  two  fit  parameters  aeff  and  io,  ae//  is  the 

89 

Distribution  Statement  A.  Approved  for  public  release;  distribution  unlimited. 

6  Of  16  88ABW-2015-5729 


American  Institute  of  Aeronautics  and  Astronautics 


r 


Figure  6.  Left:  Electric  field  propagation,  Middle:  Beer-Lambert  law  fit  to  power  density  and  SROM  results, 
Right:  Distributed  heat  source.  All  images  are  from  the  same  30%  volume  fraction  composite. 


most  important  because  it  describes  the  effective  wave  propagation  and  attenuation  through  the  depth  of  the 
composite  whereas  Iq  is  simply  an  input  to  the  model.  The  SROM  curve  in  the  middle  plot  of  Figure  6  shows 
power  distribution  if  SROM  were  used  to  calculate  ae//.  Clearly,  the  SROM  over-predicts  aeff  and  causes 
the  heat  source  to  essentially  be  a  surface  heat  source.  The  right  plot  in  Figure  6  shows  the  heat  source  as  a 
function  of  depth  for  the  FEA  simulation  and  assuming  Beer-Lambert  law  absorption.  Although  Equation  2 
fits  the  intensity  distribution  well,  using  the  extracted  aeff  to  define  the  heat  source  as  a  function  of  depth 
(Equation  3)  does  not  agree  very  well  with  the  FEA  results.  This  is  because  the  intensity  distribution  from 
the  FEA  accounts  for  absorption,  reflection,  multiple  reflections,  and  scattering  whereas  Equation  3  assumes 
the  reduction  in  intensity  is  only  due  to  absorption. 

D.  Thermal  Response  to  Optical  Energy 

Typically  the  thermal  response  to  optical  input  for  a  material  system  is  simulated  by  either  a  one  dimensional 
surface  flux  approximation  (ID  SFA)  or  a  logarithmic  heat  source.20,21  The  simulation  is  then  treated  as  a 
thermal  transport  problem  with  no  feedback  to  the  optical  solution  except  in  the  form  of  element  removal  or 
’’element  death”  in  some  FE  simulations  where  ablation  thresholds  are  reached.4  While  this  approach  may 
be  applicable  to  homogeneous  materials  or  heterogeneous  materials  with  pulsed  lasers,  it  is  not  necessarily 
valid  for  continuous  wave  (CW)  laser  interaction  with  composite  materials  having  constituents  with  varying 
optical  properties  that  lead  to  non-uniform  absorption  and  reflections.  These  varying  properties  and  non¬ 
dilute  loading  conditions  of  inclusions  on  the  order  of  the  wavelength  also  make  it  difficult  to  analytically 
predict  the  intensity  distribution  due  to  scattering  and  absorption. 

In  this  work,  determining  the  thermal  response  to  optical  input  within  a  composite  is  accomplished 
using  a  multi-fidelity  approach  within  a  single  simulation.  Two  zones  of  differing  fidelity  are  created  to 
simulate  the  optical  energy  absorption  and  the  thermal  response.  The  high  fidelity  zone  is  used  to  solve 
for  the  wave  propagation  and  absorption,  which  in  turn  becomes  the  heat  source  input  for  the  thermal 
simulation.  The  size  of  this  zone  will  depend  on  the  depth  required  to  be  considered  optically  thick,  the 
overall  part  thickness,  the  wavelength  of  interest,  and  computational  resource  constraints.  The  lower  fidelity 
zones  (“Mesh  Transition  Zone”  and  “Bulk  Material  Zone”  in  Figure  7)  are  present  to  include  the  influence  of 
convection  and  radiation  at  the  non-irradiated  (back  surface)  of  the  part  while  taking  into  account  the  actual 
part  thickness  and  different  external  thermal  conditions  on  the  irradiated  surface  (front  surface)  and  the  back 
surface.  The  geometry  created  to  encompass  both  fidelity  zones  and  the  simulation  steps  corresponding  to 
each  zone  are  shown  in  Figure  7. 

In  Figure  7,  the  spatial  heat  source  (  Qsrc )  is  determined  by  the  optical  simulation  whereas  the  effective 
specific  heat  capacity  ( Cp,eff )  and  effective  density  ( peff )  are  determined  using  SROM.  The  decision  step 
condition  of  Tr  <=  0.5%  was  determined  via  numerical  convergence  on  R,  Tr,  A1  and  ae//.  Once  the 
boundary  conditions  are  set,  Qsrc  has  been  determined,  and  the  effective  thermal  properties  have  been 
determined,  a  transient  heat  transfer  simulation  is  performed  to  determine  the  temperature  field.  All  of  the 
simulation  steps  in  this  work  are  completed  using  COMSOL.  It  should  be  noted  that  the  geometry  in  Figure 
7  makes  the  inherent  assumption  that  the  incident  optical  beam  width  is  larger  than  the  thickness  of  the 
composite.  It  should  also  be  noted  that  the  solution  most  closely  resembles  what  is  happening  directly  under 
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Figure  7.  Optical  to  thermal  algorithm  and  corresponding  geometric  configuration 


the  center  of  the  beam  and  will  not  be  as  accurate  near  the  beam  edge  due  to  the  semi-infinite  assumptions 
made  by  this  model.  In  the  areas  labeled  “Bulk  Material  Zone”  and  “Mesh  Transition  Zone”,  the  effective 
thermal  properties  are  assumed  to  be  quasi-isotropic  such  that  the  conductivity  is  isotropic  within  the  plane 
being  simulated.  For  the  purposes  of  comparing  this  multi-fidelity  approach,  a  simple  one  dimensional  heat 
transfer  model  is  constructed  which  uses  the  same  thermal  boundary  conditions  as  the  model  in  Figure  7, 
but  assumes  the  heat  source  either  comes  from  a  point  flux  at  the  surface  (ID  SFA),  or  as  a  distributed  heat 
source  as  defined  by  the  attenuation  coefficient  and  Equation  3  (ID  RMG  Input).  The  ID  RMG  Input  model 
is  a  low  fidelity  thermal  model  which  uses  the  effective  thermal  properties  and  effective  optical  properties  as 
determined  by  the  higher  fidelity  RMG  model.  The  effective  thermal  properties  for  the  ID  SFA  model  are 
found  using  the  EIM  whereas  the  effective  optical  properties  assume  all  the  incident  energy  is  absorbed  at 
the  surface. 

The  thermal  boundary  conditions  used  to  determine  the  temperature  field  of  the  opto-thermal  response 
will  greatly  impact  the  temperature  solution  inside  the  material.  In  this  work,  the  top  and  bottom  boundaries 
in  Figure  7  are  perfectly  insulated,  which  assumes  semi-infinite  conditions.  The  convection  conditions  are 
assumed  to  be  asymmetric  with  heat  transfer  coefficients  of  50  on  the  front  surface  (left  surface  in 
Figure  7),  and  10  ,^K  on  the  back  surface  (right  surface  in  Figure  7).  Surface  to  ambient  radiation  is  also 
considered  on  the  front  and  back  surfaces  because  of  the  relatively  high  (0.88)  surface  emissivity  (e)  of  the 
epoxy.  The  total  composite  thickness  used  in  this  work  is  1  mm  and  all  initial  temperatures  and  surrounding 
temperatures  are  assumed  to  be  20  °C. 

III.  Results  and  Discussion 


A.  Thermal  Conductivity 

The  ultimate  goal  of  this  work  is  to  determine  the  temperature  field  of  PMCs  after  optical  energy  absorption; 
thus  accurately  predicting  the  effective  through  thickness  thermal  conductivity  is  essential.  Although  the 
transverse  conductivity  can  be  found  in  a  number  of  ways,  the  FE  method  with  RMG  is  used  due  to 
its  versatility,  stochastic  capability,  and  lack  of  restrictions  regarding  constituent  and  interface  properties. 
Effective  thermal  conductivity  predictions  are  compared  to  experimental  results22  in  Figure  8  using  the 
following  methods:  RMG  FE,  EIM,  SROM,  HCP  unit  cell,  and  SC  unit  cell. 

The  material  system  used  to  generate  Figure  8  consists  of  pitch  based  carbon  fibers  in  an  epoxy  matrix 
where  the  fiber  conductivity  is  200  W/(m  *  K)  while  the  matrix  conductivity  is  0.3  W/(m  *  K)  to  match 
the  material  system  used  by  Thornburg  and  Pears22  and  verify  the  thermal  conductivity  model.  In  general, 
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Figure  8.  Effective  through-thickness  thermal  conductivity  predictions  compared  to  literature  values 


as  the  ratio  of  fiber  conductivity  to  matrix  conductivity  increases,  the  more  the  random  arrangement  and 
subsequent  development  of  conductive  paths  will  matter.  It  should  be  noted  that  Thornburg  and  Pears22 
observed  a  large  increase  in  porosity  in  the  two  highest  volume  fraction  composites,  which  could  account  for 
the  non-monotonic  behavior  of  the  transverse  thermal  conductivity.  As  is  shown  in  Figure  8,  the  SROM  does 
a  poor  job  of  predicting  the  effective  conductivity,  but  can  still  be  used  as  a  lower  bound  assuming  perfect 
interfaces  between  constituents.  The  other  analytic  method  included  in  Figure  8  is  the  EIM  which  agrees  very 
closely  with  the  hexagonal  unit  cell  results;  however,  both  of  these  methods  predict  lower  conductivity  than 
the  experimental  results,  RMG  FEM,  and  the  SC  unit  cell  method.  The  EIM  is  known  to  have  limitations 
on  conductivity  ratios,  which  are  exceeded  in  this  particular  case.9 

The  unit  cell  analyses  fall  victim  to  the  some  of  the  same  shortcomings  as  the  EIM  in  terms  of  neglecting 
conductive  paths  that  may  develop  due  to  the  stochastic  nature  of  the  microstructure.  Another  advantage 
of  using  the  RMG  FEM  compared  with  other  methods  is  the  ability  to  perform  stochastic  simulations  to 
find  a  range  of  properties.  Standard  error  for  each  volume  fraction  simulated  using  the  RMG  FEM  is  found 
by  generating  5  realizations  of  the  same  microstructure.  The  error  bars  included  in  Figure  8  are  of  standard 
error  and  are  quite  small  even  for  a  large  conductivity  mismatch  ratio.  Based  on  these  results  and  literature 
knowledge,  some  computational  expense  may  be  saved  by  using  an  analytic  method  such  as  the  EIM  when 
the  relative  constituent  properties,  volume  fraction,  and  interface  properties  are  within  the  allowable  range 
of  the  analytic  method.  In  general,  the  EIM  is  valid  for  nearly  all  interface  thermal  resistance  situations, 
low  volume  fraction,  and  a  small  mismatch  ratio  (<  100)  of  fiber  conductivity  to  matrix  conductivity.  The 
PMC  system  used  in  this  work  is  based  on  the  properties  of  PAN  based  carbon  fibers,  i.e.  Toray  T300™, 
and  the  epoxy  system  EPON  828™  with  thermal  conductivities  of  10.46  W/ (m  *  K )  and  0.223  W/(m  *  K) 
respectively.  The  thermal  conductivity  of  the  composite  system  as  predicted  by  various  methods  is  shown 
in  Figure  9. 

Much  like  the  comparison  in  Figure  8,  Figure  9  shows  the  RMG  method  consistently  predicting  the 
highest  thermal  conductivity  of  any  method  while  the  EIM  analytic  method  and  HCP  unit  cells  agree  very 
closely  with  each  other.  Also  similar  to  Figure  8,  Figure  9  shows  that  the  conductivity  predictions  begin 
to  diverge  at  larger  fiber  volume  fractions,  differing  by  19%  for  SC  unit  cells  and  22%  for  HCP  unit  cells 
and  EIM  at  60%  fiber  volume  fraction.  Assuming  perfect  interfaces,  the  reason  for  the  higher  conductivity 
in  the  RMG  case  is  again  because  of  the  conductive  paths  that  develop  due  to  the  random  placement. 
Unsurprisingly  the  orientation  of  the  unit  cells  has  a  negligible  effect  on  the  effective  conductivity;  however, 
unit  cell  orientation  will  prove  to  be  a  factor  when  considering  effective  optical  properties. 
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Figure  9.  Thermal  conductivity  predictions  of  the  PAN  based  carbon  fiber-epoxy  system 


B.  Optical  Properties 

Effective  optical  properties  are  useful  inputs  to  lower  fidelity  laser  interaction  models,  or  on  their  own  as 
composite  properties.  However,  the  resulting  temperature  field  can  be  found  within  the  micromechanics 
framework  without  calculating  the  effective  optical  properties  if  the  framework  is  fully  coupled.  The  most 
important  optical  properties  related  to  the  thermal  response  are  R,  Tr,  A,  and  aeff.  These  metrics  together 
describe  the  magnitude  and  location  of  the  transduction  of  optical  energy  to  thermal  energy.  In  the  interest 
of  balancing  computational  time  and  accuracy,  the  “Increase  RVE  depth”  step  shown  in  Figure  7  is  limited 
to  resizing  to  a  maximum  depth  of  150  /im.  As  shown  in  Figure  5,  150  /im  may  not  be  considered  optically 
thick  for  some  synthetic  generation  methods  at  low  volume  fractions  which  means  their  transmission  terms 
may  be  non-zero.  However,  because  of  the  method  used  to  determine  effective  attenuation  coefficient,  the 
RVE  is  not  required  to  be  optically  thick  to  determine  aef  / .  The  material  system  used  for  optical  simulation 
is  also  T300/EPON  828  with  the  electrical  and  optical  properties  summarized  in  Table  1.  The  fiber  optical 
properties  are  the  averaged  from  values  for  graphite  and  amorphous  carbon  since  PAN  based  carbon  fiber 
structure  may  be  either  of  these  depending  on  processing.23  The  electrical  properties  and  matrix  optical 
properties  are  specified  by  the  manufacturer. 

Table  1.  Electrical  and  Optical  Properties 

Property  Epoxy  Fiber 

Electrical  conductivity,  ^  10~12  58,862.5 

Real  refractive  index  1.57  2.68 

Imaginary  refractive  index  0  1.49 


The  combined  effects  of  reflection,  transmission,  absorption,  and  scattering  will  all  contribute  to  oteff. 
Considering  that  most  reflection  takes  place  near  the  surface  and  the  transmission  for  a  given  RVE  should 
be  as  close  to  zero  as  possible,  the  loss  of  intensity  as  a  function  of  depth  is  primarily  due  to  absorption 
and  scattering.  Scattering  effectively  increases  the  optical  path  length,  which  leads  to  more  absorption  if  the 
fibers  are  in  fact  absorbing.  The  predictions  for  the  effective  attenuation  coefficient  as  a  function  of  volume 
fraction  using  various  methods  is  shown  in  Figure  10. 

Figure  10  indicates  that  the  HCP  30°  unit  cell  method  most  closely  resembles  the  RMG  method  in  terms 
of  effective  attenuation  values.  Considering  Figures  5  and  10,  the  HCP  30°  unit  cell  agrees  with  the  RMG 
method  on  metrics  related  to  optical  interaction,  but  differs  when  considering  thermal  conductivity  (Figure 
9).  The  error  bars  present  on  the  RMG  results  in  Figure  10  are  the  standard  error  from  5  stochastic  iterations 
at  each  volume  fraction.  Finding  the  ae//  for  a  given  volume  fraction  requires  several  stochastic  simulations 
using  the  RMG  method  in  order  to  predict  a  value  with  reasonable  confidence.  The  SROM  curve  in  Figure 
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Figure  10.  Effective  attenuation  coefficient  (aejj)  as  a  function  of  volume  fraction  using  various  prediction 
methods. 


10  indicates  that  the  simple  analytic  solution  does  not  provide  a  reasonable  approximation  for  ae//  because 
it  basically  assumes  that  everything  is  absorbed  at  the  surface.  It  should  be  noted  that  these  assertions  only 
apply  to  opaque  fibers  in  a  transparent  matrix  where  the  wavelength  of  the  incident  light  is  smaller  than  the 
fibers.  The  polarization  dependence  of  these  claims  is  also  not  evaluated  in  this  work. 

C.  Opto-Thermal  Response 

The  temperature  held  due  to  incident  optical  energy  will  vary  as  a  function  of  several  parameters  such  as 
radiation  duration  and  intensity,  location  (depth),  thermal  boundary  conditions,  wavelength,  polarization, 
synthetic  generation  method,  and  composite  type.  Important  composite  properties  that  will  affect  the 
temperature  solution  include  the  constituent  properties  (optical  and  thermal),  relative  volume  fractions  of 
the  fibers  and  matrix,  and  the  thickness  of  the  composite  part.  Though  all  these  factors  play  a  role  in  the 
thermal  response,  only  the  effects  due  to  volume  fraction  and  generation  method  are  studied  in  this  work. 
The  fiber  volume  fraction  is  of  particular  interest  because  it  is  a  physical  property  that  impacts  both  the 
optical  properties  and  effective  thermal  properties.  The  generation  method  has  also  been  shown  to  affect 
both  the  effective  conductivity  (Figure  9)  and  the  attenuation  coefficient  (Figure  10),  however  its  effect  is 
numerical  and  should  be  mitigated.  The  motivation  to  employ  a  unit  cell  method  over  the  RMG  method 
is  to  save  computational  time.  The  thermal  response  metrics  considered  are  the  surface  temperatures  as  a 
function  of  time  and  the  average  temperature  profile  as  a  function  of  depth. 

The  temperature  as  a  function  of  time  is  directly  correlated  with  the  total  energy  absorbed  by  the 
microstructure;  therefore,  if  the  synthetic  generation  method  affects  the  overall  power  absorption  it  will 
ultimately  influence  the  thermal  response.  For  this  material  system,  Figures  5  and  10  indicate  that  mi¬ 
crostructure  generation  will  affect  temperature  more  severely  at  low  volume  fraction  (10%  -  20%)  for  SC  0°, 
HCP  0°,  and  SC  45°  methods  due  to  the  difficulty  of  creating  an  optically  thick  RVE  which  in  turn  affects 
overall  power  absorption.  The  HCP  0°  method  at  high  volume  fraction  will  also  invoke  a  different  thermal 
response  due  to  a  higher  power  absorption  coefficient  (0.879  at  60%)  than  all  the  other  2D  FE  methods 
(0.821-0.833  at  60%)  which  manifests  as  a  higher  temperature  and  faster  temperature  increase  at  both  the 
front  and  back  surfaces  of  the  composite.  The  surface  temperatures  as  a  function  of  time  are  shown  in  Figure 
11. 

The  SC  0°,  HCP  0°,  and  SC  45°  unit  cell  models  show  significantly  lower  surface  temperatures  at  low 
volume  fraction  (Figures  11A  anclllB)  due  to  their  inability  to  be  optically  thick  within  the  150  pm  depth 
constraint.  Even  at  low  volume  fractions,  the  SFA  predicts  much  higher  temperatures  at  both  the  front 
and  back  surfaces  which  is  a  direct  result  of  assuming  all  the  incident  energy  is  absorbed  at  the  surface, 
which  neglects  possible  reflections  and  multiple  reflections  due  to  the  carbon  fibers.  Effectively,  this  means 
the  SFA  assumes  more  total  energy  is  absorbed  than  the  RMG  and  unit  cell  methods.  A  similar  yet  less 
severe  effect  is  shown  in  Figures  11C  11D  for  the  HCP  0°  method  again  because  of  a  higher  overall  power 
absorption.  The  ID  RMG  Input  method  agrees  very  closely  with  the  RMG  method  for  both  high  and  low 
volume  fraction  composites  which  demonstrates  the  RMG  method’s  usefulness  as  an  input  to  lower  fidelity 
models  when  trying  to  predict  surface  temperatures.  The  position  of  the  knee  of  the  curve  in  Figures  11B 
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Figure  11.  Front  and  back  surface  temperatures  as  a  function  of  time  for  10%  volume  fraction  and  60%  volume 
fraction  using  various  prediction  methods. 


and  11D  correspond  to  a  Fourier  number  of  approximately  1  when  the  total  composite  thickness  is  the 
characteristic  length  whereas  this  distance  is  typically  the  half  width  of  a  semi-infinite  wall.  As  expected, 
the  knee  is  pushed  further  to  the  right  at  low  volume  fractions  due  to  the  lower  effective  thermal  diffusivity 
of  the  composite.  Considering  the  knee  in  the  back  surface  temperature  curves,  the  temperature  profile  plots 
shown  in  Figure  12  occur  at  2  seconds  to  ensure  enough  time  has  elapsed  to  allow  thermal  diffusion  from 
the  heat  source  in  the  RVE  near  the  front  surface  to  the  back  surface  for  all  volume  fraction  conditions. 

Figures  12A  -  12D  show  interesting  trends  in  agreement  between  the  different  methods  used.  At  very  low 
volume  fraction  (Figure  12A)  the  SC  0°,  SC  45°,  and  HCP  0°  predict  much  lower  temperatures  throughout 
the  thickness  due  to  not  being  optically  thick,  which  is  the  same  effect  seen  in  the  front  and  back  surface 
temperature  figures.  However,  the  HCP  30°,  ID  RMG  Input,  and  the  ID  SFA  mostly  agree  except  the  ID 
SFA  near  the  surface  and  the  overall  temperature  profile  shape.  The  simultaneous  agreement  at  the  back 
surface  and  disagreement  at  the  front  surface  can  be  attributed  to  the  relative  locations  of  the  heat  sources 
in  each  microstructure  and  the  asymmetric  thermal  boundary  conditions.  The  heat  source  in  the  SFA  is 
directly  at  the  surface  with  the  largest  heat  transfer  coefficient  whereas  the  heat  source  in  the  RMG  and 
HCP  30°  is  not  as  large  overall,  yet  the  sources  are  located  somewhat  below  the  surface.  Since  the  effective 
conductivities  are  approximately  equal  at  this  low  volume  fraction,  the  thermal  resistance  between  the  heat 
source  and  the  back  surface  is  approximately  equal  in  each  case.  Also  aiding  in  the  SFA  and  RMG  method 
agreeing  at  very  low  volume  fraction  is  the  power  reflection  coefficient  of  the  RMG  is  proportional  to  the 
volume  fraction  which  means  the  overall  heat  source  magnitude  is  more  comparable  between  the  RMG  and 
SFA  at  very  low  volume  fraction.  The  ID  RMG  Input  method  and  the  RMG  method  mostly  agree  at  low 
volume  fraction  except  near  the  surface.  This  difference  can  be  attributed  to  the  method  of  applying  the 
heat  source  as  shown  in  Figure  6  where  the  ID  RMG  Input  method  has  the  largest  source  at  the  surface  as 
opposed  to  the  RMG’s  largest  heat  source  which  is  generally  some  distance  below  the  surface. 

Figure  12D  shows  that  even  though  the  total  power  input  is  approximately  equal  (except  the  SFA) 
the  temperature  profile  and  temperature  gradient  varies  between  microstructure  generation  types.  The 
difference  between  the  temperature  profiles,  their  slope  in  particular,  at  high  fiber  volume  fraction  is  due 
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Figure  12.  Temperature  profiles  through  the  depth  using  various  prediction  methods 


to  the  differences  in  thermal  diffusivity,  especially  in  the  regions  with  homogenized  properties.  The  RMG 
method  always  predicts  the  highest  thermal  conductivity  which  causes  it  to  have  a  more  uniform  temperature 
through  the  thickness,  i.e.  a  higher  back  surface  temperature  and  a  lower  front  surface  temperature  than 
methods  with  comparable  total  power  input.  The  ID  RMG  Input  model  agrees  very  well  with  the  RMG 
model  at  high  volume  fraction  because  the  high  density  of  fibers  forces  the  largest  part  of  the  distributed 
heat  source  closer  to  the  surface,  which  is  what  the  ID  RMG  Input  model  assumes.  However,  unlike  the  ID 
SFA  model,  the  ID  RMG  Input  model  also  accounts  for  total  power  reflection  and  uses  the  higher  effective 
thermal  conductivity  values  calculated  by  the  RMG  model.  In  comparing  the  ID  SFA  model  to  the  ID 
RMG  Input  model,  it  is  clear  that  accounting  for  reflection,  a  distributed  heat  source,  and  correct  effective 
thermal  conductivity  is  necessary  to  predict  the  thermal  response  with  a  low  fidelity  model. 

At  moderate  volume  fractions  (  30%  and  40%  in  Figures  12B  and  12C)  the  various  prediction  methods 
agree  well  with  each  other  aside  from  the  ID  SFA  predicting  somewhat  higher  temperature  due  to  a  larger 
total  heat  source.  The  agreement  at  these  volume  fractions  represents  a  balance  of  two  competing  numerical 
mechanisms,  optical  absorption  and  effective  thermal  conductivity.  The  trade-off  between  thermal  conduc¬ 
tivity  agreement  and  attenuation  coefficient  agreement  as  a  function  of  volume  fraction  is  shown  in  Figure 
13. 

Figure  13  shows  better  thermal  conductivity  agreement  at  lower  volume  fractions  for  various  generation 
types  whereas  they  begin  to  diverge  at  larger  volume  fractions.  Meanwhile,  it  is  difficult  to  achieve  an 
optically  thick  RVE  with  good  values  for  attenuation  coefficient  at  low  volume  fraction  with  unit  cell  methods 
but  easy  to  achieve  this  condition  at  high  volume  fraction.  Thus,  the  various  generation  methods  agree  on 
thermal  conductivity  at  low  volume  fraction  and  agree  on  optical  absorption  at  high  volume  fraction.  This 
trade-off  is  why  the  various  generation  methods  do  not  agree  at  volume  fraction  extremes  while  agreeing 
much  better  at  moderate  volume  fractions. 
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Figure  13.  Relative  agreement  with  RMG  effective  properties  at  various  volume  fractions  for  various  generation 
methods. 


IV.  Conclusion 

We  have  developed  a  stochastic,  micromechanics  framework  in  Matlab  and  COMSOL  to  find  the  thermal 
response  of  a  PMC  due  to  incident  optical  energy.  The  micromechanics  framework  developed  in  this  work 
finds  the  effective  thermal  properties,  effective  optical  properties,  heat  source  and  thermal  response  to  optical 
energy  within  the  microstructure  in  a  stochastic  manner.  We  have  shown  that  commonly  employed  analytic 
models  cannot  always  sufficiently  predict  effective  thermal  properties  and  are  not  valid  for  optical  property 
prediction.  This  work  has  also  demonstrated  the  ability  to  solve  the  full  wave  solution  of  Maxwell’s  equations 
inside  the  random  and  ordered  microstructure  of  a  composite  which  accounts  for  diffraction  and  interference 
effects  that  are  neglected  by  geometric  optics  and  surface  flux  approximations. 

Effective  thermal  properties  can  generally  be  found  using  analytic  expressions  or  unit  cell  analysis  in 
many  composite  systems,  but  it  was  shown  that  this  is  not  the  case  all  composite  systems  and  accurate 
prediction  of  the  effective  conductivity  is  vital  to  predicting  the  opto-thermal  response.  The  applicability 
of  this  model  was  demonstrated  in  a  composite  with  high  thermal  conductivity  mismatch  ratio  where  the 
series  resistance  between  the  fibers  (a  function  of  fiber  spacing  and  matrix  conductivity)  plays  a  large  role. 
Although  unit  cell  analysis  can  typically  be  used  to  find  effective  thermal  properties,  it  does  not  perform 
well  when  finding  effective  optical  properties  due  to  the  fundamental  differences  between  thermal  diffusion 
and  wave  propagation. 

The  periodic  nature  of  the  unit  cell  method  allows  clear  propagation  paths  at  lower  volume  fractions 
within  the  microstructure  that  are  highly  unlikely  to  arise  in  real  composites.  These  paths  diminish  the 
contribution  of  the  fiber  properties  to  the  effective  optical  properties  when  the  matrix  is  transparent.  Should 
the  optical  properties  of  the  matrix  and  the  fibers  be  reversed,  the  effect  of  the  periodic  arrangement  of 
the  fibers  may  not  have  as  much  impact  on  the  results.  Thus,  for  the  case  of  highly  absorbing  fibers  in  a 
transparent  matrix,  a  random  microstructure  must  be  used  as  opposed  to  a  periodic  microstructure  in  order 
to  find  effective  optical  properties,  especially  the  total  power  input  and  the  effective  attenuation  coefficient. 

The  thermal  response  to  optical  energy  takes  into  account  both  the  effective  thermal  properties  and 
the  optical  absorption  of  the  composite.  We  have  assessed  the  applicability  of  using  a  unit  cell  analysis,  ID 
thermal  analysis,  and  a  random  RVE  to  find  the  thermal  response  of  a  composite  due  to  absorption  of  optical 
energy.  At  low  volume  fractions,  the  unit  cell  models  predict  lower  temperatures  than  the  RMG  model  due  to 
their  higher  transmission  (smaller  attenuation  coefficient)  through  the  RVE,  thus  reducing  the  heat  source. 
However,  at  higher  fiber  volume  fraction  the  unit  cell  methods  predict  higher  temperatures  on  the  front 
surface  and  lower  temperatures  at  the  back  surface  when  compared  to  the  RMG  method.  This  is  due  to 
the  effective  attenuation  coefficients  and  total  power  input  being  approximately  the  same  but  the  effective 
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conductivity  is  higher  when  using  the  RMG  method  which  allows  more  efficient  dissipation  of  the  thermal 
energy  through  the  thickness  of  the  composite.  Two  ID  transient  thermal  models  were  also  evaluated  for 
their  ability  to  predict  the  temperature  profile  through  the  thickness.  It  was  shown  that  the  SFA  assuming 
that  all  energy  is  absorbed  at  the  surface  leads  to  much  higher  temperature  predictions  than  all  the  other 
methods.  This  is  primarily  due  to  neglecting  reflections  and  assuming  the  heat  source  is  concentrated  at  the 
surface  instead  of  being  applied  as  a  distributed  source. 

The  agreement  between  the  ID  RMG  Input  model  and  the  RMG  model  demonstrates  the  usefulness  of  this 
micromechanics  framework  and  its  ability  to  find  effective  optical  and  thermal  properties  which  can  subse¬ 
quently  be  used  as  inputs  to  lower  fidelity  models.  The  developed  framework  predicts  the  microscale  thermal 
response  of  composites  under  laser  irradiation  for  various  composite  properties  which  allows  optimization  of 
processing  steps  such  as  laser  curing  and  thermoplastic  consolidation  where  the  desired  temperature  response 
is  within  a  certain  range.  The  full  wave  solution  within  the  microstructure  includes  the  effects  of  fiber  ar¬ 
rangement  and  other  microscale  inclusions  that  can  significantly  affect  laser  absorption  and  subsequent  heat 
transfer.  These  microscale  effects,  which  can  manifest  as  a  macro  scale  temperature  response,  are  neglected 
by  lower  fidelity  models  and  are  not  fully  captured  by  geometric  optics  approximations. 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 


AFM 

atomic  force  micrograph 

AFRL 

Air  Force  Research  Laboratory 

ALD 

atomic  layer  deposition 

BMI 

bismaleinmide 

BMPM 

bismaleimidodiphenyl  methane 

FEA 

Finite  element  analysis 

FEM 

finite  element  models 

ICMSE 

Integrated  computational  materials  science  and  engineering 

MD 

molecular  dynamics 

OMC 

organic  matrix  composite 

SThM 

scanning  thermal  microscopy 

WPAFB 

Wright-Patterson  Air  Force  Base 
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